Molecular Dynamics Investigation of Solution Structure between NaCl

ARTICLE pubs.acs.org/JPCC

Molecular Dynamics Investigation of Solution Structure between NaCl and Quartz Crystals Melanie B. Webb,† Stephen H. Garofalini,‡ and George W. Scherer*,§ †

Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States Interfacial Molecular Science Laboratory, Department of Materials Science and Engineering, Rutgers University, Piscataway, New Jersey 08855, United States § Department of Civil and Environmental Engineering, Princeton University, Princeton, New Jersey 08544, United States ‡

ABSTRACT: Molecular dynamics simulations with a dissociative interatomic potential were used to examine the interaction forces between quartz and sodium chloride crystals in the presence of pure water and a supersaturated (7 m NaCl) solution at crystal separations between 7 and 24.5 Å. The force between the crystals was found to be attractive for all separations at both ion concentrations. The structuring of interfacial water molecules, leading to ordering of the dipole moments at both interfaces, is primarily responsible for the attraction. Modifications of the interfacial water structure by the decreasing crystal separation are more prevalent at the quartz/liquid interface than at the sodium chloride/water interface.

’ INTRODUCTION Forces between solids immersed in a liquid are important to many research areas including the stability of colloids,1 flotation and separation of minerals,2 pushing of particles by growing crystals,3
5 and crystallization pressure.6,7 When such forces are repulsive and occur between a crystal growing in confinement and the confining body (such as can occur within a porous material), damaging stresses can arise. Such is the case with destructive crystallization pressure by salt or ice, which is responsible for sculpting mountains and damaging monuments.8 However, in some systems, the net interaction is attractive. It has been reported that quartz and sodium chloride are such a system, where the attraction has been attributed to the modification of bulk solution structure through the orientation of the water dipoles at the surface.2 This is important for mineral flotation and also for understanding the potential for damage to quartzitic sandstone monuments exposed to marine salts. In this study, we use molecular dynamics to confirm the nature of the interaction between halite and quartz in the presence of pure water or a slightly supersaturated solution. We also examine the relative influence of crystal separation on the structuring of the solution near the two crystal surfaces. We first describe the potential model, the development of the systems of study, and the computational simulation method. The following section shows the effects of concentration and crystal separation on the structuring of the confined solutions and pressure felt by the confined sodium chloride crystal.

both two-body and three-body terms. The form of the pair potential used for all interactions is Upair ¼ Uqq þ Uqd qd þ Uqd q þ Uqqd þ Urep þ Udisp

ð1Þ where

r 2011 American Chemical Society

ð2Þ

!
 j rij rij qid qd erf erfc Uqd qd ðrij Þ ¼ rij 2ξij β

ð3Þ

!
 qid qj rij rij erf pffiffiffi erfc Uqd q ðrij Þ ¼ rij β 2ξij

ð4Þ

!
 j rij rij qi qd erf pffiffiffi erfc Uqqd ðrij Þ ¼ rij β 2ξij

ð5Þ

Urep ðrij Þ ¼

2Aijrep ξijr rij

rij erfc 2ξijr

! ð6Þ

ij

Udisp ðrij Þ ¼

’ COMPUTATIONAL PROCEDURE Potential Model. A fully atomistic potential developed by Mahadevan and Garofalini for dissociative water9 and water
silica studies,10 which allows for bonds to be broken or formed between all atoms, was used in this work. This potential contains


 qi qj rij erfc Uqq ðrij Þ ¼ rij β


C6 rij 6

ð7Þ

Received: May 3, 2011 Revised: August 24, 2011 Published: August 27, 2011 19724

dx.doi.org/10.1021/jp204093c | J. Phys. Chem. C 2011, 115, 19724–19732

The Journal of Physical Chemistry C

ARTICLE

Table 1. Potential Parameters (a) Parameters of the Two-body Potential species

Arep (J)

ξ (Å)

ξr (Å)

C6 (J Å6)

Na
Na11 Na
Cl11 Cl
Cl11 Na
O11 Na
H11 Cl
O11 Cl
H11 O
H9 O
O9 H
H9 Si
O10 Si
Si10 Si
H10 Si
Cl Si
Na

1.998  10
19 2.111  10
17 7.063  10
19 2.685  10
17 5.787  10
19 3.626  10
17 1.831  10
18 2.283  10
16 4.250  10
17

24 24 24 24 24 24 24 24 24 24 24 24 24 24 24

1.004 0.6434 0.997 0.5375 0.6171 0.6952 0.5745 f(T,P)a 0.610

1.698  10
19 9.925  10
19 4.170  10
18 1.015  10
18

0.373 0.640 0.350 0.3868 0.5225

7.000  10
18

2.670  10
16 7.000  10
17 5.000  10
16 7.177  10
16 2.839  10
17

4.821  10
18

4.226  10
18

3.800  10
18

(b) Charges on Species q/e +1.0
1.0
0.904 +0.452 +1.808

species Na11 Cl11 O9 H9 Si10

qd/e
0.25 +0.25 +0.226
0.113
0.452

(c) Three-Body Parameters species H
O
H9 Na
O
Na11 Na
O
H11

O
Si
O10 Si
O
Si10 Si
O
H10

Na
O
Si

λjik (ergs) 30  10
11 1  10
11 5  10
11 ij = O
Na ik = O
H 15  10
11 1  10
11 5  10
11 ij = O
Si ik = O
H 1  10
11 ij = O
Si ik = O
Na

r0ij (Å) 1.6 2.8

γij (Å) 1.3 2.0

2.8 1.5 3.0 2.8

2.0 1.2 2.8 2.0

2.8 1.5

2.0 1.2

2.8 2.8

2.0 2.0

θjik (o) 100 109.5 109.5

109.5 109.5 109.5

109.5

(d) Columns B(:0)
B(:2) of matrix Bmn in ξOH (T,P), defined in footnote a r 0.655726502 3.403472  10
4
4.057853  10
8 1.657262  10
12


1.04442689  10
2
3.986929  10
6 4.677537  10
10
1.838785  10
14

8.31892416  10
5 1.742261  10
8
2.007873  10
12 7.549619  10
17

(e) Columns B(:3)
B(:5) of matrix Bmn in ξOH (T,P), defined in footnote a r
3.07929142  10
7
3.364186  10
11 3.800411  10
15
1.355453  10
19 a

5.44770929  10
10 2.419996  10
14
2.672717  10
18 8.939302  10
23


3.73609493  10
13 0 0 0

m=3,n=5 m n ξOH r (T,P) = ∑m=0,n=0 BmnP T for 263 e T(K) e 373 and 1 e P (atm) e 6000.

19725

dx.doi.org/10.1021/jp204093c |J. Phys. Chem. C 2011, 115, 19724–19732

The Journal of Physical Chemistry C

ARTICLE

As in previous studies using this potential,9
11 qdi and qi in eqs 3
5 are related through the relationship qdi = qi/4. Equations 2
5 are the Coulombic terms and each contains a complementary error function component that serves as a dampening function on these long-range interactions. This is part of the Wolf summation12 that was used instead of the Ewald sum in this work. The Coulombic energy contribution, with the appropriate choice of the dampening parameter, β, and cut off distance, Rc, is Eele

1 N ¼ 2 i¼1

∑ ∑ j6¼ i

qi qj erfcðrij =βÞ
rij

lim

rij sfRc

qi qj erfcðrij =βÞ rij

!!

rij < Rc





 N erfcðRc =βÞ 1 þ qi 2 2Rc βπ1=2 i ¼ 1

∑

ð8Þ

A cutoff radius of 1 nm and dampening parameter of 4.46 Å were used in the Wolf summation for the long-range Coulombic interactions. A three-body potential of the following form was applied to the interaction potential: U3 ðri , rj , rk Þ ¼ v3 ðrij , rik , θjik Þ

ð9Þ

where

! γik ðcosðθjik Þ
cosðθ0jik ÞÞ2 v3 ðrij , rik , θjik Þ ¼ λjik exp þ rij
rij0 rik
rik0 γij

ð10Þ when rik < r0ik and rij < r0ij, and otherwise υ3(rij0 rik0 θjik) = 0. The three-body interactions only exist for certain ijk triplets present in the system and are identically zero for all others. This term serves as an energy penalty when three atoms form an angle that deviates from the ideal angle θ0jik. The parameters for the pair and three-body interactions used in this work, except for Si
Cl and Si
Na interactions, were determined previously9
11 and are shown in parts a
e of Table 1. All MD simulations were performed using a fifth-order Nordsiek-Gear predictor-corrector algorithm with a time step of 0.1 fs, unless otherwise noted. The simulations were performed in the NVE (constant number, volume, and energy), NVT (constant number, volume, and temperature), or NPT (constant number, pressure, and temperature) ensembles. Parameter Determination. The Si
Cl and Si
Na interaction parameters were determined to replicate the structures of a SiCl4 molecule and sodium silicate glass, when used in conjunction with the other previously determined parameters. A sodium silicate glass was created by randomly placing 200 Si atoms, 450 O atoms, and 100 Na atoms into a box, such that the density was 2.393 g/cm3 at 300 K. This system was put through a meltquench process at temperatures of 300, 6000, 5000, 4000, 3000, 2000, 1000, 300 K under NVT conditions for the first 20 ps followed by 80 ps under NVE conditions for each temperature using a time step of 1 fs. Using the parameters in Table 1, this system yielded a structure with pair distribution function peaks centered at 1.63 Å, 3.15 Å, and 2.64 Å for Si
O, Si
Si, and O
O respectively and broad peaks centered at 3.32 Å for Si
Na and 2.31 Å for O
Na. These peak positions are similar to those found in previous simulation and experimental studies, as summarized in Newell et al.13

Figure 1. Initial configuration of a 7 m quartz-solution
sodium chloride system. Solution in front of and behind the NaCl crystal in the x direction has been removed so that the crystal can be seen clearly. Oxygen atoms are red, silicon are yellow, hydrogen are white, sodium are blue, and chlorine are teal. The dashed area represents the gap area of interest and the shaded region indicates the immobile quartz atoms. Snapshot was created using VMD.15

The Si
Cl interaction parameters were determined using a SiCl4 molecule at 1 and 200 K. Four Cl atoms were placed tetrahedrally around a central Si atom and run for 5 ps under NVT conditions followed by 15 ps under NVE conditions without the Wolf summation because this system is an isolated molecule. The resultant configuration was cooled to 1 K with 5 ps of NVT simulation followed by another 25 ps under NVE conditions. The last 10 ps of the 1 K simulation yielded a Si
Cl distance 2.020 Å, whereas the 200 K simulation yielded 2.028 Å, which is comparable to the 2.019 Å found in electron diffraction studies.14 System of Study. To study the effects of concentration and gap distance on the structure and behavior of the solution between quartz and sodium chloride crystal surfaces, a sodium chloride crystal was placed at varying distances from the quartz surfaces to produce slit-type gaps filled by either pure water (0 m) or a supersaturated 7 m NaCl solution. An illustration, created with VMD,15 of a starting configuration can be found in Figure 1. This was accomplished by first relaxing an 11  14  7 unit cell (9702 atoms) of α-quartz, oriented such that the z axis is normal to the (001) plane, under NPT conditions at 298 K and 1 atm for 50 ps with a time step of 1 fs. An equilibrium volume was reached within the first 2.5 ps of the simulation and varied with a standard deviation of 0.05% of the average value subsequently. The resultant crystal was 56.2  61.96  39.3 Å and was shifted slightly in the z direction, so that one surface (A) had only Si atoms exposed while the other (B) had only O atoms exposed. The relaxed quartz crystal was then prehydroxylated by placing 2902 water molecules randomly above the xy plane of the quartz with periodic boundaries to permit reactions on both the A and B surfaces of the quartz. A relaxation simulation at constant volume, consisting of 2.5 ps in NVT at 300 K followed by 7.5 ps in NVE, was carried out to allow water molecules to begin to react with the surface and allow the pressure caused by the random initial placement of the water molecules to moderate before the constant volume restriction was removed from the z 19726

dx.doi.org/10.1021/jp204093c |J. Phys. Chem. C 2011, 115, 19724–19732

The Journal of Physical Chemistry C direction. During the NVE portion, the average temperature of the system varied within 4 K of the desired 298 K. Reaction was allowed to continue for 0.15 ns under NPT conditions, where only the box size in the z direction responds to changes in system pressure. During the simulation, water reacts with the quartz surfaces, leaving SiOH and SiOH2 groups on the surfaces and hydroniums, hydroxyls, and unreacted water in the solution. The solution was then removed from the system, leaving only the mostly hydroxylated quartz to be used for the next stage. In addition to the quartz, a 6  7  5 unit cell (1680 atoms) sodium chloride was relaxed under NPT conditions at 298 K and 1 atm for 50 ps with a 1 fs time step. The crystal reached an equilibrium volume after 31 ps and varied with a standard deviation of 0.02% from that point forward. The resultant NaCl crystal had (001) planes on the surfaces and dimensions of 33.84  39.48  28.20 Å. The sodium chloride crystal was placed in a system of size 56.2  61.96  66.2 Å centered in the xy plane and at z positions ranging from 2.5 to 20 Å from the base of the box. Pure water or ∼7 m NaCl solution was randomly placed in the unoccupied volume of the box, so as to achieve a density at 25 °C of 0.997 g/cm3 for water16 and 1.22 g/cm3 for 7 m NaCl (extrapolated from aqueous sodium chloride solution densities17). This amounted to 6428 water molecules for pure water, whereas 5585 water molecules and 704 NaCl ion pairs were added to the total volume for the 7 m NaCl solution. Special attention was given to make sure that the concentration of NaCl in solution directly beneath the NaCl crystal was as close to 7 m as possible because this was the gap area examined in this study. The gap was filled with Na+ and Cl
ions and water molecules at the appropriate density before the remaining ions and water were placed into the rest of the surrounding area. The sodium chloride crystal and solution system was then stacked on top of the relaxed and prehydroxylated quartz with approximately 1.8 Å of space between the hydroxylated surfaces and the solution surrounding the NaCl crystal. The total gap distance between the quartz surface A and the NaCl crystal surface is defined as the z distance between the highest Si atom on surface A and the lowest Na or Cl atom in the NaCl crystal, and was held constant during our simulations. In addition to the water molecules and Na and Cl ions, 5 hydronium ions, and 170 hydroxyl ions that were removed from the hydroxylated quartz system were reintroduced into the total system to ensure overall charge neutrality. Simulation Details. The design of our system allows solution to move into and out of the gap space, replenishing or diminishing the concentration or density, as the system requires. As we wished to examine the effects of both concentration and gap height on the structure of the solution, we fixed the distance between the highest Si atom on Quartz Surface A and the lowest Na or Cl in the sodium chloride crystal during the course of the simulation. After each move in the simulation, the highest Si atom in Quartz Surface A was found and the atoms within the sodium chloride crystal were moved en masse, such that it was exactly the original gap distance away. This was further facilitated by not allowing the atoms within the sodium chloride crystal to move under the forces of the other system atoms, such that the lowest sodium chloride atom in the crystal is the same throughout the simulation. This choice was made to prevent any dissolution or rearrangement of the atoms within the crystal and to preserve the crystal shape and size so that the x and y dimensions of the gap space were precisely known. Furthermore,

ARTICLE

3080 central atoms located more than 12 Å from either quartz surface in the quartz crystal were immobile. Both of these immobile quartz atoms and the immobile sodium chloride crystal atoms were allowed to react to changes in system pressure during constant pressure simulations. Each simulation consisted of a short constant-volume relaxation simulation, which allowed the randomly placed atoms and molecules to adjust into more natural positions without causing the overall system size to react wildly to the rapidly changing forces, followed by three longer constant pressure simulations. In all cases, periodic boundary conditions were used in all directions. The relaxation portion of the simulation consisted of 5 ps under NVT conditions at a temperature of 298 K followed by 20 ps under NVE conditions. During the three constant-pressure portions of the simulations, only the z dimension of the simulation box could respond to the changing system pressure. This allowed the x and y dimensions of the box to remain the same but still produce a system in which the density of the solution could change. Therefore, the gap height between surface B of the quartz crystal and the top surface of the NaCl crystal varied throughout the length of the simulation. The first section of the constant pressure simulations consisted of 0.1 ns of simulation time under NPT conditions of 298 K and 1 atm. This was followed by two 0.2 ns simulations under the same conditions, which gave a total simulation time of 0.5 ns. During these portions of the simulation, the z force on the NaCl crystal caused by the solution and quartz crystal underneath it was recorded. All analysis was performed on the final 0.2 ns of these simulations.

’ RESULTS AND DISCUSSION Structure. Structuring of the solution within the gap between the crystal surfaces was examined using atom density profiles in the z direction from the final 0.2 ns of the simulations. Parts a
f of Figure 2 contains a representative set (hb ≈ 24.5, 16.5, 8.5 Å) of these profiles for both the 0 and 7 m systems. In each of the profiles, the quartz crystal is located on the left side of the plot (smaller z position values), whereas the NaCl crystal is located on the right side, at higher z. In all of the systems, the Si layer of the quartz crystal ends at z ≈ 40 Å. Additionally, a silanol layer is observed starting ∼2 Å within the quartz crystal and extending ∼2 Å past the Si surface of the original quartz crystal. There are also H+ attached to oxygens that bridge between two Si atoms within the quartz. Two layers of these bridging OH are observed around the first two oxygen peaks within the quartz, ∼5 Å from the surface. These were likely formed during the reaction of water with the surface Si atoms forming a silanol and an H+, which penetrated several layers within the crystal and attached to a bridging oxygen. Protonated bridging oxygens in a vitreous silica and water system have been observed in MD18 and molecular orbital19,20 calculations and have been previously discussed.18 Furthermore, a small layer of water molecules is attached to the Si on the quartz surface (Si
O distance