Supporting Information for âMolecular Simulation of Carbon Dioxide ...
Supporting Information for “Molecular Simulation of Carbon Dioxide, Brine, and Clay Mineral Interactions and Determination of Contact Angles” Craig M. Tenney∗ and Randall T. Cygan Sandia National Laboratories, Albuquerque, New Mexico, USA E-mail: [email protected] Phone: 505-284-4466. Fax: 505-844-7786 This Supporting Information document contains the following sections: • Simulation Details: details describing the setup and analysis of the systems studied in this work • Force Fields: parameters used for the calculation of energies and forces in the molecular simulations described in this work • Example LAMMPS Scripts: examples of the LAMMPS scripts used to set up and conduct the molecular simulations described in this work
S1
Figure S1: Crystal structure of a kaolinite unit cell. Oxygen atoms are depicted by red vertices, hydrogen by white, aluminum by pink, and silicon by yellow.
Simulation Details Figure S1 shows the structure for a kaolinite unit cell. Each silicon is connected to three other silicon atoms by bridging oxygens forming the siloxane surface. This arrangement of tetrahedral silicons gives rise to hexagonal cavities on the surface. In the gibbsite sheet, each aluminum is coordinated by two oxygens of the siloxane sheet and shares four hydroxyls with neighboring aluminum atoms. Three of these hydroxyls are oriented toward the external surface forming the gibbsite surface while the fourth one is oriented inward in the direction of the siloxane cavity. This layered structure extends to edges of the crystal where exposed dangling silicon and aluminum atoms are terminated by hydroxyls at low pH conditions. The mineral phase used in each simulation was composed of three layers of kaolinite. Initial atom positions within the kaolinite slab were assigned by replicating a unit cell derived from single crystal synchrotron data (1 ). The experimental unit cell shape was modified to allow use of orthogonal simulation cells. Built from 38 x 22 x 3 unit cells (85272 atoms), the kaolinite slab measures 195.9 Å and 196.7 Å in the transverse (x and y) directions and 19.4 Å thick. Simulation cell boundaries in the x and y directions are periodic, fixed, and sized to match the dimensions of the kaolinite slab, resulting in an infinite slab with no edges. In the ∗
To whom correspondence should be addressed
S2
z direction simulation cells are non-periodic and sized to contain the kaolinite slab and the fluid adjacent to one side of the slab. The fluid region adjacent to the slab is approximately 140 Å thick. Kaolinite was modeled using the fully-flexible ClayFF force field (2 ). ClayFF uses an empirically derived set of interaction parameters to accurately describe the potential energy between atoms in the clay structure. ClayFF has been used extensively to successfully simulate many oxide, hydroxide, and hydrated systems including bulk and interfacial structures (2 –6 ). Metal-oxygen interactions are described by a 12-6 Lennard-Jones (LJ) term and a Coulombic function with partial charges derived by quantum chemistry methods. The only explicit bonded interactions are those within hydroxyl groups (bond stretch), which are based on the SPC water model (7 ). To improve computational efficiency and also prevent unphysical oscillation of the three-layer kaolinite slab, atoms in the two kaolinite layers furthest from the fluid interface were held fixed during the simulations. The water and dissolved ion models used in this work were chosen for consistency with ClayFF. ClayFF uses the three-point simple point charge (SPC) water model (7 ), combined with harmonic bond stretching and angle bending terms based on the intramolecular parameters from Teleman et al.(8 ). Interatomic potentials for water that incorporate bond flexibility have demonstrated improved predictions of properties across a wide range of state points compared to rigid models (9 , 10 ). Having one LJ center at the oxygen atom and partial charges centered directly on each of the three atoms, the SPC water model is relatively simple and has been used to study properties of bulk water and aqueous systems (7 , 8 , 11 , 12 ). ClayFF uses force field parameters from the literature to represent the aqueous species Ca2+ (13 ), Na+ (14 ), and Cl– (14 ). Supercritical CO2 was also modeled using a fully flexible force field (15 ). In this threepoint force field, each atom is represented by a LJ center and partial charge to account for non-bonded energies. Bond stretch and angle bend terms account for intramolecular energies. The force field uses the LJ parameters of Zhu et al. (16 ), which were refined from
S3
the original model of Harris and Yung (17 ). The harmonic intramolecular bond stretch and angle bend parameters were optimized to reproduce the vibrational spectra of CO2 (15 ). All simulations were conducted using the open-source LAMMPS (version: 25 Jul 2012) molecular dynamics software package (18 , 19 ). Long range electrostatic interactions are handled in a pseudo 2D manner by a routine in LAMMPS that imposes a large vacuum space between periodic images in the z direction and corrects for dipole interactions between images. Temperature was maintained via a Nose-Hoover thermostat implemented in the standard LAMMPS code. Pressure was controlled by a virtual piston (implemented in custom LAMMPS code) attached to a repulsive wall located adjacent to the fluid surface opposite and parallel to the kaolinite slab. Initial system configurations were all constructed using Packmol (20 ), an open-source simulation setup tool. Total simulation times of 10 ns to 14 ns were used to allow droplets to equilibrate with the surrounding bulk fluid before collecting data for analysis. For a typical system in this study (∼500k atoms), 1 ns of simulation time required approximately 12 hours of wall clock time, using 64 dual socket/quad core (2.93 GHz, Nehalem X5570) compute nodes, i.e. 512 cores. System trajectory snapshots were recorded every 10 ps. These trajectories were postprocessed to visualize molecular configurations and to calculate average local densities of CO2 , water, and brine species. A system was assumed to be sufficiently equilibrated when the number and distribution of molecules within the droplet stopped changing with time. Contact angles were calculated assuming spherical geometry at the fluid–fluid interface. Contact angles θ for wetting droplets were calculated from measurements of droplet height h and droplet base radius b according to
θ = arcsin
2hb + h2
b2
(1)
Contact angles for non-wetting droplets were calculated as the average of the two independent
S4
values calculated according to b r
(2)
h−r r
(3)
θ = π − arcsin and θ = π − arccos
where r is the droplet radius. The goals of this project did not justify the development of routines for automated calculation of contact angles from MD trajectories, so droplet dimensions h, b, and r were estimated manually from plots of molecule density.
Force Fields FORCEFIELD CO2-Cygan
! SETTING: !! vdw_scaling: van der Waals interaction scaling coeffs for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! vdw_cutoff: !!! style (str, (’cut’, ’shift’, ’correct’)): cutoff style !!! cutoff (float, angstrom): cutoff distance !! vdw_switch: !!! style (str, (’charmm’,)): switch style !!! switch_inner (float, angstrom): inner switch distance !!! switch_outer (float, angstrom): outer switch distance !! vdw_mix: !!! style (str, (’LJ’,)): interaction style
S5
!!! mixrule (str, (’LB’, ’geometric’)): mixing rule !! charge_scaling: charge interaction scaling coefficients for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! charge_cutoff: !!! style (str, (’cut’, ’ewald’)): cutoff style !!! cutoff (float, angstrom): cutoff distance # typically for use with CLAYFF SETTING
! 5
vdw_scaling 0.0000 0.0000 1.0000 vdw_cutoff cut 12.0 # ? vdw_mix LJ LB charge_scaling 0.0000 0.0000 1.0000 charge_cutoff ewald 12.0 # ?
! ATOM: !! name (str): atom name (required) !! mass (float, amu): atomic mass !! element (str): element symbol !! epsilon (float, Ken): LJ well depth !! sigma (float, angstrom): LJ zero-energy distance (rmin/2**(1/6)) !! charge (float, e): partial atomic charge ATOM name charge epsilon sigma element mass ! 2 c
0.6512
28.144 2.800
C
12.011 # CO2 carbon ! 1
o
-0.3256
80.378 3.028
O
15.999 # CO2 oxygen ! 2
S6
! BOND: !! harmonic: K * (R - R0)**2 !!! K (float, Ken/angstrom^2): force constant !!! R0 (float, angstrom): equilibrium bond length !! fixed: R = R0 !!! R0 (float, angstrom): equilibrium bond length BOND
! 1
c
o harmonic
1015458
1.162
! 1
1.160
! ANGLE: !! harmonic: K * (A - A0)**2 !!! K (float, Ken/radian^2): force constant !!! A0 (float, degrees): equilibrium angle !! fixed: A = A0 !!! A0 (float, degrees): equilibrium angle ANGLE o
! 1 c
o harmonic
54351
180.00
MOLECULE CO2
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate
S7
! 1
180.000
ATOM name type x y z ! 3 C
c
0.0000
0.0000
0.0000
! 1
O1
o
-1.1600
0.0000
0.0000
! 2
O2
o
1.1600
0.0000
0.0000
! 3
! BOND: BOND
! 2
C
O1
! 1
1.160
C
O2
! 2
1.160
! ANGLE: ANGLE
! 1
O1
C
O2
! 1
180.000
# converted from OpenMD CLAYFF.frc FORCEFIELD H2O-clayff
! SETTING: !! vdw_scaling: van der Waals interaction scaling coeffs for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! vdw_cutoff: !!! style (str, (’cut’, ’shift’, ’correct’)): cutoff style !!! cutoff (float, angstrom): cutoff distance !! vdw_switch: !!! style (str, (’charmm’,)): switch style
S8
!!! switch_inner (float, angstrom): inner switch distance !!! switch_outer (float, angstrom): outer switch distance !! vdw_mix: !!! style (str, (’LJ’,)): interaction style !!! mixrule (str, (’LB’, ’geometric’)): mixing rule !! charge_scaling: charge interaction scaling coefficients for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! charge_cutoff: !!! style (str, (’cut’, ’ewald’)): cutoff style !!! cutoff (float, angstrom): cutoff distance SETTING
! 5
vdw_scaling 0.0000 0.0000 1.0000 vdw_cutoff cut 12.0 # ? vdw_mix LJ LB charge_scaling 0.0000 0.0000 1.0000 charge_cutoff ewald 12.0 # ?
! ATOM: !! name (str): atom name (required) !! mass (float, amu): atomic mass !! element (str): element symbol !! epsilon (float, Ken): LJ well depth !! sigma (float, angstrom): LJ zero-energy distance (rmin/2**(1/6)) !! charge (float, e): partial atomic charge ATOM name charge epsilon sigma element mass ! 2
S9
h*
0.41
0 0.000
H
o*
-0.82
78.25256 3.166
O
1.008 # water hydrogen ! 1 15.999 # water oxygen ! 2
! BOND: !! harmonic: K * (R - R0)**2 !!! K (float, Ken/angstrom^2): force constant !!! R0 (float, angstrom): equilibrium bond length !! fixed: R = R0 !!! R0 (float, angstrom): equilibrium bond length BOND h*
! 1 o* harmonic
279038
1.000
! 1
0.969
! ANGLE: !! harmonic: K * (A - A0)**2 !!! K (float, Ken/radian^2): force constant !!! A0 (float, degrees): equilibrium angle !! fixed: A = A0 !!! A0 (float, degrees): equilibrium angle ANGLE
! 1
h*
o*
h* harmonic
23048
109.47
MOLECULE H2O
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge
S10
! 1
103.996
!! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate ATOM name type x y z ! 3 O
o*
0.0000
-0.0643
0.0000
! 1
H1
h*
-0.7633
0.5321
0.0000
! 2
H2
h*
0.7633
0.5321
0.0000
! 3
! BOND: BOND
! 2
H1
O
! 1
0.969
H2
O
! 2
0.969
! ANGLE: ANGLE
! 1
H1
O
H2
! 1
103.996
# converted from OpenMD CLAYFF.frc # This is the forcefield file for the Clay Force Field (CLAYFF) # Details can be found in the following article: # "Molecular Models of Hydroxide, Oxyhydroxid, and Clay Phases and # the Development of a General Force Field" by Randall T. Cygan, # Jian-Jie Liang, and Andrey G. Kalinichev, J. Phys. Chem. B 108, # pp. 1255-1266 (2004). FORCEFIELD clayff
! SETTING:
S11
!! vdw_scaling: van der Waals interaction scaling coeffs for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! vdw_cutoff: !!! style (str, (’cut’, ’shift’, ’correct’)): cutoff style !!! cutoff (float, angstrom): cutoff distance !! vdw_switch: !!! style (str, (’charmm’,)): switch style !!! switch_inner (float, angstrom): inner switch distance !!! switch_outer (float, angstrom): outer switch distance !! vdw_mix: !!! style (str, (’LJ’,)): interaction style !!! mixrule (str, (’LB’, ’geometric’)): mixing rule !! charge_scaling: charge interaction scaling coefficients for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! charge_cutoff: !!! style (str, (’cut’, ’ewald’)): cutoff style !!! cutoff (float, angstrom): cutoff distance SETTING
! 5
vdw_scaling 0.0000 0.0000 1.0000 vdw_cutoff cut 12.0 # ? vdw_mix LJ LB charge_scaling 0.0000 0.0000 1.0000 charge_cutoff ewald 12.0 # ?
S12
! ATOM: !! name (str): atom name (required) !! mass (float, amu): atomic mass !! element (str): element symbol !! epsilon (float, Ken): LJ well depth !! sigma (float, angstrom): LJ zero-energy distance (rmin/2**(1/6)) !! charge (float, e): partial atomic charge ATOM name charge epsilon sigma element mass ! 24 h*
0.41
0 0.000
H
1.008 # water hydrogen ! 1
ho
0.425
0 0.000
H
1.008 # hydroxyl hydrogen ! 2
o*
-0.82
78.25256 3.166
O
15.999 # water oxygen ! 3
oh
-0.95
78.25256 3.166
O
15.999 # hydroxyl oxygen ! 4
ob
-1.05
78.25256 3.166
O
15.999 # bridging oxygen ! 5
obos
-1.1808
78.25256 3.166
O
15.999 # bridging oxygen with octahedral subst
obts
-1.16875
78.25256 3.166
O
15.999 # bridging oxygen with tetrahedral subs
obss
-1.2996
78.25256 3.166
O
15.999 # bridging oxygen with double substitut
ohs
-1.0808
78.25256 3.166
O
15.999 # hydroxyl oxygen with substitution ! 9
st
2.1 0.0009267943 3.302 Si
28.085 # tetrahedral silicon ! 10
ao
1.575 0.0006696284 4.271 Al
26.982 # octahedral aluminum ! 11
at
1.575 0.0009267943 3.302 Al
26.982 # tetrahedral aluminum ! 12
mgo
1.36 0.0004547008 5.264 Mg
24.305 # octahedral magnesium ! 13
mgh
1.05 0.0004547008 5.264 Mg
24.305 # hydroxide magnesium ! 14
cao
1.36
0.002532785 5.567 Ca
40.078 # octahedral calcium ! 15
cah
1.05
0.002532785 5.562 Ca
40.078 # hydroxide calcium ! 16
feo
1.575
0.004547008 4.906 Fe
55.845 # octahedral iron ! 17
lio
0.525
0.004547008 4.210 Li
S13
6.941 # octahedral lithium ! 18
Na
1
65.5126 2.350 Na
K
1
50.35557 3.334
Cs
1
50.35557 3.831 Cs
Ca
2
50.35557 2.872 Ca
Ba
2
23.66711 3.817 Ba
Cl
-1
50.40593 4.400 Cl
K
22.990 # aqueous sodium ion ! 19 39.098 # aqueous potassium ion ! 20 132.905 # aqueous cesium ion ! 21 40.078 # aqueous calcium ion ! 22 137.327 # aqueous barium ion ! 23 35.453 # aqueous chloride ion ! 24
! BOND: !! harmonic: K * (R - R0)**2 !!! K (float, Ken/angstrom^2): force constant !!! R0 (float, angstrom): equilibrium bond length !! fixed: R = R0 !!! R0 (float, angstrom): equilibrium bond length BOND
! 3
h*
o* harmonic
279038
1.000
! 1
ho
oh harmonic
279038
1.000
! 2
ho
ohs harmonic
279038
1.000
! 3
# ANGLE # h*
o*
h* harmonic 23048 109.47
# ao
oh
ho harmonic 15107 109.47
# ao
ohs ho harmonic 15107 109.47
# at
oh
# at
ohs ho harmonic 15107 109.47
# mgh oh
ho harmonic 15107 109.47
ho harmonic 15107 109.47
# mgh ohs ho harmonic 15107 109.47 # cah oh
ho harmonic 15107 109.47
S14
0.794
# cah ohs ho harmonic 15107 109.47 # feo oh
ho harmonic 15107 109.47
# feo ohs ho harmonic 15107 109.47 # lio oh
ho harmonic 15107 109.47
# lio ohs ho harmonic 15107 109.47 MOLECULE kaolinite
! SETTING: !! unitcell: unit cell dimensions (for crystalline materials only) !!! a (float, angstrom): unit cell edge length !!! b (float, angstrom): unit cell edge length !!! c (float, angstrom): unit cell edge length !!! alpha (float, degree): unit cell angle !!! beta (float, degree): unit cell angle !!! gamma (float, degree): unit cell angle !!! space_group (str, P1): space group SETTING
! 1
unitcell
5.154
8.942
7.401
91.69 104.61
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate ATOM name type x y z ! 34
S15
89.82 P1
Al1
ao
0.6650
4.3300
3.4040
! 1
Al2
ao
3.2140
2.8550
3.3960
! 2
Si1
st
4.9750
3.0050
0.6610
! 3
Si2
st
2.4620
1.4720
0.6710
! 4
O1
ob
-0.3210
3.0970
2.2630
! 5
O2
ob
0.0550
5.8590
2.2660
! 6
O3
ob
0.0140
4.4710
0.0000
! 7
O4
ob
1.0450
2.0680
0.1750
! 8
O5
ob
1.0710
6.8310
0.0020
! 9
O6
oh
-0.3200
8.5930
2.3290
! 10
O7
oh
3.8190
1.3530
4.3260
! 11
O8
oh
-0.9220
4.1030
4.3240
! 12
O9
oh
-0.9230
7.5290
4.3520
! 13
H1
ho
0.0860
0.2420
2.4870
! 14
H2
ho
-1.0150
1.4610
5.0180
! 15
H3
ho
-1.1230
4.1950
5.0680
! 16
H4
ho
-1.1110
6.9610
4.9970
! 17
Al3
ao
3.2560
8.8010
3.4040
! 18
Al4
ao
0.6510
7.3260
3.3960
! 19
Si3
st
2.4120
7.4760
0.6610
! 20
Si4
st
-0.1010
5.9430
0.6710
! 21
O10
ob
2.2700
7.5680
2.2630
! 22
O11
ob
2.6180
1.3880
2.2660
! 23
O12
ob
2.5770
0.0000
0.0000
! 24
O13
ob
3.6360
6.5390
0.1750
! 25
O14
ob
3.6340
2.3600
0.0020
! 26
O15
oh
2.2430
4.1220
2.3290
! 27
S16
O16
oh
1.2560
5.8240
4.3260
! 28
O17
oh
1.6690
8.5740
4.3240
! 29
O18
oh
1.6400
3.0580
4.3520
! 30
H5
ho
2.6770
4.7130
2.4870
! 31
H6
ho
1.5760
5.9320
5.0180
! 32
H7
ho
1.4690
8.6660
5.0680
! 33
H8
ho
1.4520
2.4900
4.9970
! 34
! BOND: BOND
! 8
H1
O6
! 1
0.750
H2
O7
! 2
0.770
H3
O8
! 3
0.776
H4
O9
! 4
0.880
H5
O15
! 5
0.750
H6
O16
! 6
0.770
H7
O17
! 7
0.776
H8
O18
! 8
0.880
# for use with CLAYFF MOLECULE Na
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate
S17
!! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate ATOM name type x y z ! 1 Na
Na
0.0000
0.0000
0.0000
! 1
# for use with CLAYFF MOLECULE Ca
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate ATOM name type x y z ! 1 Ca
Ca
0.0000
0.0000
0.0000
! 1
# for use with CLAYFF MOLECULE Cl
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate
S18
!! z (float, angstrom): atom coordinate ATOM name type x y z ! 1 Cl
Cl
0.0000
0.0000
0.0000
! 1
Example LAMMPS Scripts System Setup # read lammps data file, write a restart file
# these variables must be set from the command line: variable
outname index init
variable
infile index init.data
# optional # initial temperature (do nothing if < 0) variable
beg_temp index -1
# these variables have generally conservative defaults: # kspace: 1 = pppm, 2 = ewald, 3 = ewald/n, other = off variable
kspace index 1
# restart != 0: write_restart at end of script variable
restart index 1
# ___end command line variable section___
# set up simulation # the following information is saved to restart files S19
units
real
atom_style
full
boundary
p p f
pair_style
lj/cut/coul/long 12
pair_modify
tail no
pair_modify mix arithmetic special_bonds lj 0.000000 0.000000 1.000000 special_bonds coul 0.000000 0.000000 1.000000 bond_style harmonic angle_style harmonic
read_data
${infile}
# end information that is saved to restart files
if "${kspace} == 1" then "kspace_style pppm 0.0001" if "${kspace} == 2" then "kspace_style ewald 0.0001" if "${kspace} == 3" then "kspace_style ewald/n 0.0001" kspace_modify slab 3.0
if "${beg_temp} >= 0" then & "velocity all create ${beg_temp} 24601 mom yes rot yes dist gaussian"
if "${restart} != 0" then & "write_restart
${outname}.restart.*"
S20
System Relaxation # NVE, NVT, NPT, or NPH ensemble
# these variables must be set variable outname index init.relax variable infile index init.restart.0 # ensemble/integrator: 1 = NVE, 2 = NVT, 3 = NPT, 4 = NPH, 5 = NVT/relax variable ensemble index 5 variable numsteps index 1000000
# damping terms are in multiples of ${timestep} # thermostat settings (if applicable) variable beg_temp index 1 variable end_temp index 320 variable tdamp index 10 # barostat settings (if applicable) variable beg_press index 200 variable
end_press index ${beg_press}
variable
pdamp index 1000
# ptensor = iso | aniso | tri | x | y | z | xy | xz | yz variable
ptensor index iso
# output options: # thermo_freq: output thermo snapshots at specified interval variable
thermo_freq index 1000
# dump_freq: dump trajectory snapshots at specified interval variable dump_freq index 1000 S21
# dump_detail: 1 = positions, 2 = +velocities, 3 = +forces variable
dump_detail index 1
# press_freq > 0: output pressure tensor data at specified interval variable
press_freq index 0
# other options: # kspace: 1 = pppm, 2 = ewald, 3 = ewald/n, other = off # (e.g set kspace=0 for non-periodic systems) variable
kspace index 1
# timestep: aka dt, e.g. 1 fs variable
timestep index 1
# processor allocation (i.e. "processors ${procx} ${procy} ${procz}") variable
procx index *
variable
procy index *
variable
procz index *
# restart_freq: save restart file every N steps (0 = only at end, = 0: reset current step to new_step variable
new_step index -1
# density > 0: rescale volume to result in specified density (avoid NPT, NPH) # note for ’real’ units: g/cm^3 * 0.6022 = g/mol/A^3 variable
density index 0
# ___end command line variable section___
# set up simulation
S22
processors ${procx} ${procy} ${procz} read_restart
${infile}
if "${kspace} == 1" then "kspace_style pppm 0.0001" if "${kspace} == 2" then "kspace_style ewald 0.0001" if "${kspace} == 3" then "kspace_style ewald/n 0.0001" neighbor
2.0 bin
neigh_modify
delay 5 check yes
timestep
${timestep}
kspace_modify slab 3.0
# change this to hold some atoms rigid #group MOBILE union all region MOBILE block INF INF INF INF 148 162 side out units box group MOBILE region MOBILE group FIXED subtract all MOBILE neigh_modify exclude group FIXED FIXED delete_bonds FIXED multi remove
# initial wall position should match previous final wall position variable zwall index "zlo 10" fix Zwall MOBILE wall/harmonic ${zwall} 50 0 1 units box extpress ${beg_press} 1e-8 variable FtoP equal 68568.411/lx/ly
# atm / (kcal/mol / Angstrom)
variable Pwall equal f_Zwall[1]*${FtoP}
# reset current simulation timestep if "${new_step} >= 0" then "reset_timestep ${new_step}"
S23
# pick desired ensemble/integrator variable _tdamp equal ${tdamp}*${timestep} variable _pdamp equal ${pdamp}*${timestep} if "${ensemble} == 1" then & "fix NVE MOBILE nve" & elif "${ensemble} == 2" & "fix TFIX MOBILE nvt temp ${beg_temp} ${end_temp} ${_tdamp}" & elif "${ensemble} == 3" &
"fix TFIX MOBILE npt temp ${beg_temp} ${end_temp} ${_tdamp} ${ptensor} ${beg_press} ${ elif "${ensemble} == 4" &
"fix TFIX MOBILE nph iso ${beg_press} ${end_press} ${_pdamp} nreset 1000 dilate partia elif "${ensemble} == 5" & "fix NVE MOBILE nve/limit 0.1" & "fix TFIX MOBILE langevin ${beg_temp} ${end_temp} ${_tdamp} 24601 tally yes" & else & "print ’Invalid ensemble/integrator!!!’" "exit"
# set up thermo data output if "${ensemble} > 1" then & "thermo_style custom step temp epair emol etotal press vol f_TFIX v_Pwall f_Zwall[2]" & else & "thermo_style custom step temp epair emol etotal press vol" thermo ${thermo_freq} thermo_modify flush yes
# set up dump/trajectory output -- potentially huge file!
S24
if "${dump_detail} == 3" then &
"dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu vx vy vz fx if "${dump_detail} == 2" then &
"dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu vx vy vz ma if "${dump_detail} == 1" then & "dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu mass q" dump_modify DUMP flush yes
# pressure tensor if "${press_freq} 0" then "restart ${restart_freq} ${outname}.restart" run
${numsteps}
if "${restart_freq} >= 0" then "write_restart ${outname}.restart.*"
NPT Simulation # NVE, NVT, NPT, or NPH ensemble
# naming conventions variable basename index t_330-p_20 variable suffix index 0 variable
outname index ${basename}.${suffix}
variable
infile index ${basename}.restart.*
# ensemble/integrator: 1 = NVE, 2 = NVT, 3 = NPT, 4 = NPH, 5 = NVT/relax variable ensemble index 2 variable numsteps index 7000000
# damping terms are in multiples of ${timestep} # thermostat settings (if applicable) variable beg_temp index 330 variable
end_temp index ${beg_temp}
variable
tdamp index 100
# barostat settings (if applicable) variable beg_press index 197.38465 variable
end_press index ${beg_press} S26
variable
pdamp index 1000
# ptensor = iso | aniso | tri | x | y | z | xy | xz | yz variable
ptensor index iso
# initial wall position should match previous final wall position variable zwall index "zlo 11.561567"
# output options: # thermo_freq: output thermo snapshots at specified interval variable
thermo_freq index 10000
# dump_freq: dump trajectory snapshots at specified interval variable
dump_freq index ${thermo_freq}
# dump_detail: 1 = positions, 2 = +velocities, 3 = +forces variable
dump_detail index 1
# press_freq > 0: output pressure tensor data at specified interval variable
press_freq index 0
# other options: # kspace: 1 = pppm, 2 = ewald, 3 = ewald/n, other = off # (e.g set kspace=0 for non-periodic systems) variable
kspace index 1
# timestep: aka dt, e.g. 1 fs variable
timestep index 1
# processor allocation (i.e. "processors ${procx} ${procy} ${procz}") variable
procx index *
variable
procy index *
variable
procz index *
# restart_freq: save restart file every N steps (0 = only at end, = 0: reset current step to new_step variable
new_step index -1
# density > 0: rescale volume to result in specified density (avoid NPT, NPH) # note for ’real’ units: g/cm^3 * 0.6022 = g/mol/A^3 variable
density index 0
# ___end command line variable section___
# set up simulation processors ${procx} ${procy} ${procz} read_restart
${infile}
if "${kspace} == 1" then "kspace_style pppm 0.0001" if "${kspace} == 2" then "kspace_style ewald 0.0001" if "${kspace} == 3" then "kspace_style ewald/n 0.0001" neighbor
2.0 bin
neigh_modify
delay 5 check yes
timestep
${timestep}
#change_box all z final -10 172 units box kspace_modify slab 3.0 diff ad
# change this to hold some atoms rigid #group MOBILE union all region MOBILE block INF INF INF INF 148 162 side out units box group MOBILE region MOBILE
S28
group FIXED subtract all MOBILE neigh_modify exclude group FIXED FIXED delete_bonds FIXED multi remove
fix Zwall MOBILE wall/harmonic ${zwall} 50 0 1 units box extpress ${beg_press} 1e-7 variable FtoP equal 68568.411/lx/ly
# atm / (kcal/mol / Angstrom)
variable Pwall equal f_Zwall[1]*${FtoP}
# reset current simulation timestep if "${new_step} >= 0" then "reset_timestep ${new_step}"
# pick desired ensemble/integrator variable _tdamp equal ${tdamp}*${timestep} variable _pdamp equal ${pdamp}*${timestep} if "${ensemble} == 1" then & "fix NVE MOBILE nve" & "compute TFIX_temp MOBILE temp" & elif "${ensemble} == 2" & "fix TFIX MOBILE nvt temp ${beg_temp} ${end_temp} ${_tdamp}" & elif "${ensemble} == 3" &
"fix TFIX MOBILE npt temp ${beg_temp} ${end_temp} ${_tdamp} ${ptensor} ${beg_press} ${ elif "${ensemble} == 4" &
"fix TFIX MOBILE nph iso ${beg_press} ${end_press} ${_pdamp} nreset 1000 dilate partia elif "${ensemble} == 5" & "fix NVE MOBILE nve/limit 0.1" & "fix TFIX MOBILE langevin ${beg_temp} ${end_temp} ${_tdamp} 24601 tally yes" & "compute TFIX_temp MOBILE temp" &
S29
else & "print ’Invalid ensemble/integrator!!!’" "exit"
# set up thermo data output if "${ensemble} > 1" then & "thermo_style custom step temp epair emol etotal press vol f_TFIX v_Pwall f_Zwall[2]" & "thermo_modify temp TFIX_temp" & else & "thermo_style custom step temp epair emol etotal press vol" thermo ${thermo_freq} thermo_modify flush yes
# set up dump/trajectory output -- potentially huge file! if "${dump_detail} == 3" then &
"dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu vx vy vz fx if "${dump_detail} == 2" then &
"dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu vx vy vz ma if "${dump_detail} == 1" then & "dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu mass q" dump_modify DUMP flush yes
# pressure tensor if "${press_freq} 0" then "restart ${restart_freq} ${basename}.restart.*" run
${numsteps}
if "${restart_freq} >= 0" then "write_restart ${basename}.restart.*"
Acknowledgement This material is based upon work supported as part of the Center for Frontiers of Subsurface Energy Security, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0001114. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DEAC04-94AL85000.
S31
References (1) Neder, R. B.; Burghammer, M.; Grasl, T.; Schulz, H.; Bram, A.; Fiedler, S. Refinement of the Kaolinite Structure from Single-Crystal Synchrotron Data. Clays Clay Miner. 1999, 47, 487–494. (2) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. Molecular Models of Hydroxide, Oxyhydroxide, and Clay Phases and the Development of a General Force Field. J. Phys. Chem. B 2004, 108, 1255–1266. (3) Cygan, R. T.; Greathouse, J. A.; Heinz, H.; Kalinichev, A. G. Molecular Models and Simulations of Layered Materials. J. Mater. Chem. 2009, 19, 2470–2481. (4) Thyveetil, M.-A.; Coveney, P. V.; Greenwell, H. C.; Suter, J. L. Computer Simulation Study of the Structural Stability and Materials Properties of DNA-Intercalated Layered Double Hydroxides. J. Am. Chem. Soc. 2008, 130, 4742–4756. (5) Wang, J. W.; Kalinichev, A. G.; Kirkpatrick, R. J.; Cygan, R. T. Structure, Energetics, and Dynamics of Water Adsorbed on the Muscovite (001) Surface: A Molecular Dynamics Simulation. J. Phys. Chem. B 2005, 109, 15893–15905. (6) Greathouse, J. A.; Cygan, R. T. Molecular Dynamics Simulation of Uranyl(VI) Adsorption Equilibria onto an External Montmorillonite Surface. Phys. Chem. Chem. Phys. 2005, 7, 3580–3586. (7) Berendsen, H.; Postma, J.; Van Gunsteren, W.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermolecular forces 1981, 11, 331–342. (8) Teleman, O.; Jonsson, B.; Engstrom, S. A Molecular-Dynamics Simulation of a Water Model with Intramolecular Degrees of Freedom. Mol. Phys. 1987, 60, 193–203. (9) Raabe, G.; Sadus, R. J. Molecular Dynamics Simulation of the Dielectric Constant of Water: the Effect of Bond Flexibility. J. Chem. Phys. 2011, 134, 6. S32
(10) Raabe, G.; Sadus, R. J. Molecular Dynamics Simulation of the Effect of Bond Flexibility on the Transport Properties of Water. J. Chem. Phys. 2012, 137, 8. (11) Wallqvist, A.; Teleman, O. Properties of Flexible Water Models. Mol. Phys. 1991, 74, 515–533. (12) Smith, D. E.; Haymet, A. D. J. Simulations of Aqueous-Solutions - the Role of Flexibility and the Treatment of Long-Range Forces. Fluid Phase Equilib. 1993, 88, 79–87. (13) Koneshan, S.; Rasaiah, J. C.; Lynden-Bell, R. M.; Lee, S. H. Solvent Structure, Dynamics, and Ion Mobility in Aqueous Solutions at 25 Degrees C. J. Phys. Chem. B 1998, 102, 4193–4204. (14) Smith, D. E.; Dang, L. X. Computer-Simulations of Nacl Association in Polarizable Water. J. Chem. Phys. 1994, 100, 3757–3766. (15) Cygan, R. T.; Romanov, V. N.; Myshakin, E. M. Molecular Simulation of Carbon Dioxide Capture by Montmorillonite Using an Accurate and Flexible Force Field. J. Phys. Chem. C 2012, 116, 13079–13091. (16) Zhu, A.; Zhang, X.; Liu, Q.; Zhang, Q. A Fully Flexible Potential Model for Carbon Dioxide. Chin. J. Chem. Eng. 2009, 17, 268–272. (17) Harris, J. G.; Yung, K. H. Carbon Dioxides Liquid-Vapor Coexistence Curve and Critical Properties As Predicted by a Simple Molecular-Model. J. Phys. Chem. 1995, 99, 12021–12024. (18) Plimpton, S. J. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. (19) Plimpton, S. J. LAMMPS Molecular Dynamics Simulator. http://lammps.sandia. gov.
S33
(20) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157–2164.
S34
S1
Figure S1: Crystal structure of a kaolinite unit cell. Oxygen atoms are depicted by red vertices, hydrogen by white, aluminum by pink, and silicon by yellow.
Simulation Details Figure S1 shows the structure for a kaolinite unit cell. Each silicon is connected to three other silicon atoms by bridging oxygens forming the siloxane surface. This arrangement of tetrahedral silicons gives rise to hexagonal cavities on the surface. In the gibbsite sheet, each aluminum is coordinated by two oxygens of the siloxane sheet and shares four hydroxyls with neighboring aluminum atoms. Three of these hydroxyls are oriented toward the external surface forming the gibbsite surface while the fourth one is oriented inward in the direction of the siloxane cavity. This layered structure extends to edges of the crystal where exposed dangling silicon and aluminum atoms are terminated by hydroxyls at low pH conditions. The mineral phase used in each simulation was composed of three layers of kaolinite. Initial atom positions within the kaolinite slab were assigned by replicating a unit cell derived from single crystal synchrotron data (1 ). The experimental unit cell shape was modified to allow use of orthogonal simulation cells. Built from 38 x 22 x 3 unit cells (85272 atoms), the kaolinite slab measures 195.9 Å and 196.7 Å in the transverse (x and y) directions and 19.4 Å thick. Simulation cell boundaries in the x and y directions are periodic, fixed, and sized to match the dimensions of the kaolinite slab, resulting in an infinite slab with no edges. In the ∗
To whom correspondence should be addressed
S2
z direction simulation cells are non-periodic and sized to contain the kaolinite slab and the fluid adjacent to one side of the slab. The fluid region adjacent to the slab is approximately 140 Å thick. Kaolinite was modeled using the fully-flexible ClayFF force field (2 ). ClayFF uses an empirically derived set of interaction parameters to accurately describe the potential energy between atoms in the clay structure. ClayFF has been used extensively to successfully simulate many oxide, hydroxide, and hydrated systems including bulk and interfacial structures (2 –6 ). Metal-oxygen interactions are described by a 12-6 Lennard-Jones (LJ) term and a Coulombic function with partial charges derived by quantum chemistry methods. The only explicit bonded interactions are those within hydroxyl groups (bond stretch), which are based on the SPC water model (7 ). To improve computational efficiency and also prevent unphysical oscillation of the three-layer kaolinite slab, atoms in the two kaolinite layers furthest from the fluid interface were held fixed during the simulations. The water and dissolved ion models used in this work were chosen for consistency with ClayFF. ClayFF uses the three-point simple point charge (SPC) water model (7 ), combined with harmonic bond stretching and angle bending terms based on the intramolecular parameters from Teleman et al.(8 ). Interatomic potentials for water that incorporate bond flexibility have demonstrated improved predictions of properties across a wide range of state points compared to rigid models (9 , 10 ). Having one LJ center at the oxygen atom and partial charges centered directly on each of the three atoms, the SPC water model is relatively simple and has been used to study properties of bulk water and aqueous systems (7 , 8 , 11 , 12 ). ClayFF uses force field parameters from the literature to represent the aqueous species Ca2+ (13 ), Na+ (14 ), and Cl– (14 ). Supercritical CO2 was also modeled using a fully flexible force field (15 ). In this threepoint force field, each atom is represented by a LJ center and partial charge to account for non-bonded energies. Bond stretch and angle bend terms account for intramolecular energies. The force field uses the LJ parameters of Zhu et al. (16 ), which were refined from
S3
the original model of Harris and Yung (17 ). The harmonic intramolecular bond stretch and angle bend parameters were optimized to reproduce the vibrational spectra of CO2 (15 ). All simulations were conducted using the open-source LAMMPS (version: 25 Jul 2012) molecular dynamics software package (18 , 19 ). Long range electrostatic interactions are handled in a pseudo 2D manner by a routine in LAMMPS that imposes a large vacuum space between periodic images in the z direction and corrects for dipole interactions between images. Temperature was maintained via a Nose-Hoover thermostat implemented in the standard LAMMPS code. Pressure was controlled by a virtual piston (implemented in custom LAMMPS code) attached to a repulsive wall located adjacent to the fluid surface opposite and parallel to the kaolinite slab. Initial system configurations were all constructed using Packmol (20 ), an open-source simulation setup tool. Total simulation times of 10 ns to 14 ns were used to allow droplets to equilibrate with the surrounding bulk fluid before collecting data for analysis. For a typical system in this study (∼500k atoms), 1 ns of simulation time required approximately 12 hours of wall clock time, using 64 dual socket/quad core (2.93 GHz, Nehalem X5570) compute nodes, i.e. 512 cores. System trajectory snapshots were recorded every 10 ps. These trajectories were postprocessed to visualize molecular configurations and to calculate average local densities of CO2 , water, and brine species. A system was assumed to be sufficiently equilibrated when the number and distribution of molecules within the droplet stopped changing with time. Contact angles were calculated assuming spherical geometry at the fluid–fluid interface. Contact angles θ for wetting droplets were calculated from measurements of droplet height h and droplet base radius b according to
θ = arcsin
2hb + h2
b2
(1)
Contact angles for non-wetting droplets were calculated as the average of the two independent
S4
values calculated according to b r
(2)
h−r r
(3)
θ = π − arcsin and θ = π − arccos
where r is the droplet radius. The goals of this project did not justify the development of routines for automated calculation of contact angles from MD trajectories, so droplet dimensions h, b, and r were estimated manually from plots of molecule density.
Force Fields FORCEFIELD CO2-Cygan
! SETTING: !! vdw_scaling: van der Waals interaction scaling coeffs for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! vdw_cutoff: !!! style (str, (’cut’, ’shift’, ’correct’)): cutoff style !!! cutoff (float, angstrom): cutoff distance !! vdw_switch: !!! style (str, (’charmm’,)): switch style !!! switch_inner (float, angstrom): inner switch distance !!! switch_outer (float, angstrom): outer switch distance !! vdw_mix: !!! style (str, (’LJ’,)): interaction style
S5
!!! mixrule (str, (’LB’, ’geometric’)): mixing rule !! charge_scaling: charge interaction scaling coefficients for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! charge_cutoff: !!! style (str, (’cut’, ’ewald’)): cutoff style !!! cutoff (float, angstrom): cutoff distance # typically for use with CLAYFF SETTING
! 5
vdw_scaling 0.0000 0.0000 1.0000 vdw_cutoff cut 12.0 # ? vdw_mix LJ LB charge_scaling 0.0000 0.0000 1.0000 charge_cutoff ewald 12.0 # ?
! ATOM: !! name (str): atom name (required) !! mass (float, amu): atomic mass !! element (str): element symbol !! epsilon (float, Ken): LJ well depth !! sigma (float, angstrom): LJ zero-energy distance (rmin/2**(1/6)) !! charge (float, e): partial atomic charge ATOM name charge epsilon sigma element mass ! 2 c
0.6512
28.144 2.800
C
12.011 # CO2 carbon ! 1
o
-0.3256
80.378 3.028
O
15.999 # CO2 oxygen ! 2
S6
! BOND: !! harmonic: K * (R - R0)**2 !!! K (float, Ken/angstrom^2): force constant !!! R0 (float, angstrom): equilibrium bond length !! fixed: R = R0 !!! R0 (float, angstrom): equilibrium bond length BOND
! 1
c
o harmonic
1015458
1.162
! 1
1.160
! ANGLE: !! harmonic: K * (A - A0)**2 !!! K (float, Ken/radian^2): force constant !!! A0 (float, degrees): equilibrium angle !! fixed: A = A0 !!! A0 (float, degrees): equilibrium angle ANGLE o
! 1 c
o harmonic
54351
180.00
MOLECULE CO2
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate
S7
! 1
180.000
ATOM name type x y z ! 3 C
c
0.0000
0.0000
0.0000
! 1
O1
o
-1.1600
0.0000
0.0000
! 2
O2
o
1.1600
0.0000
0.0000
! 3
! BOND: BOND
! 2
C
O1
! 1
1.160
C
O2
! 2
1.160
! ANGLE: ANGLE
! 1
O1
C
O2
! 1
180.000
# converted from OpenMD CLAYFF.frc FORCEFIELD H2O-clayff
! SETTING: !! vdw_scaling: van der Waals interaction scaling coeffs for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! vdw_cutoff: !!! style (str, (’cut’, ’shift’, ’correct’)): cutoff style !!! cutoff (float, angstrom): cutoff distance !! vdw_switch: !!! style (str, (’charmm’,)): switch style
S8
!!! switch_inner (float, angstrom): inner switch distance !!! switch_outer (float, angstrom): outer switch distance !! vdw_mix: !!! style (str, (’LJ’,)): interaction style !!! mixrule (str, (’LB’, ’geometric’)): mixing rule !! charge_scaling: charge interaction scaling coefficients for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! charge_cutoff: !!! style (str, (’cut’, ’ewald’)): cutoff style !!! cutoff (float, angstrom): cutoff distance SETTING
! 5
vdw_scaling 0.0000 0.0000 1.0000 vdw_cutoff cut 12.0 # ? vdw_mix LJ LB charge_scaling 0.0000 0.0000 1.0000 charge_cutoff ewald 12.0 # ?
! ATOM: !! name (str): atom name (required) !! mass (float, amu): atomic mass !! element (str): element symbol !! epsilon (float, Ken): LJ well depth !! sigma (float, angstrom): LJ zero-energy distance (rmin/2**(1/6)) !! charge (float, e): partial atomic charge ATOM name charge epsilon sigma element mass ! 2
S9
h*
0.41
0 0.000
H
o*
-0.82
78.25256 3.166
O
1.008 # water hydrogen ! 1 15.999 # water oxygen ! 2
! BOND: !! harmonic: K * (R - R0)**2 !!! K (float, Ken/angstrom^2): force constant !!! R0 (float, angstrom): equilibrium bond length !! fixed: R = R0 !!! R0 (float, angstrom): equilibrium bond length BOND h*
! 1 o* harmonic
279038
1.000
! 1
0.969
! ANGLE: !! harmonic: K * (A - A0)**2 !!! K (float, Ken/radian^2): force constant !!! A0 (float, degrees): equilibrium angle !! fixed: A = A0 !!! A0 (float, degrees): equilibrium angle ANGLE
! 1
h*
o*
h* harmonic
23048
109.47
MOLECULE H2O
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge
S10
! 1
103.996
!! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate ATOM name type x y z ! 3 O
o*
0.0000
-0.0643
0.0000
! 1
H1
h*
-0.7633
0.5321
0.0000
! 2
H2
h*
0.7633
0.5321
0.0000
! 3
! BOND: BOND
! 2
H1
O
! 1
0.969
H2
O
! 2
0.969
! ANGLE: ANGLE
! 1
H1
O
H2
! 1
103.996
# converted from OpenMD CLAYFF.frc # This is the forcefield file for the Clay Force Field (CLAYFF) # Details can be found in the following article: # "Molecular Models of Hydroxide, Oxyhydroxid, and Clay Phases and # the Development of a General Force Field" by Randall T. Cygan, # Jian-Jie Liang, and Andrey G. Kalinichev, J. Phys. Chem. B 108, # pp. 1255-1266 (2004). FORCEFIELD clayff
! SETTING:
S11
!! vdw_scaling: van der Waals interaction scaling coeffs for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! vdw_cutoff: !!! style (str, (’cut’, ’shift’, ’correct’)): cutoff style !!! cutoff (float, angstrom): cutoff distance !! vdw_switch: !!! style (str, (’charmm’,)): switch style !!! switch_inner (float, angstrom): inner switch distance !!! switch_outer (float, angstrom): outer switch distance !! vdw_mix: !!! style (str, (’LJ’,)): interaction style !!! mixrule (str, (’LB’, ’geometric’)): mixing rule !! charge_scaling: charge interaction scaling coefficients for bonded atoms !!! scale_12 (float): scale factor for 1-2 interactions !!! scale_13 (float): scale factor for 1-3 interactions !!! scale_14 (float): scale factor for 1-4 interactions !! charge_cutoff: !!! style (str, (’cut’, ’ewald’)): cutoff style !!! cutoff (float, angstrom): cutoff distance SETTING
! 5
vdw_scaling 0.0000 0.0000 1.0000 vdw_cutoff cut 12.0 # ? vdw_mix LJ LB charge_scaling 0.0000 0.0000 1.0000 charge_cutoff ewald 12.0 # ?
S12
! ATOM: !! name (str): atom name (required) !! mass (float, amu): atomic mass !! element (str): element symbol !! epsilon (float, Ken): LJ well depth !! sigma (float, angstrom): LJ zero-energy distance (rmin/2**(1/6)) !! charge (float, e): partial atomic charge ATOM name charge epsilon sigma element mass ! 24 h*
0.41
0 0.000
H
1.008 # water hydrogen ! 1
ho
0.425
0 0.000
H
1.008 # hydroxyl hydrogen ! 2
o*
-0.82
78.25256 3.166
O
15.999 # water oxygen ! 3
oh
-0.95
78.25256 3.166
O
15.999 # hydroxyl oxygen ! 4
ob
-1.05
78.25256 3.166
O
15.999 # bridging oxygen ! 5
obos
-1.1808
78.25256 3.166
O
15.999 # bridging oxygen with octahedral subst
obts
-1.16875
78.25256 3.166
O
15.999 # bridging oxygen with tetrahedral subs
obss
-1.2996
78.25256 3.166
O
15.999 # bridging oxygen with double substitut
ohs
-1.0808
78.25256 3.166
O
15.999 # hydroxyl oxygen with substitution ! 9
st
2.1 0.0009267943 3.302 Si
28.085 # tetrahedral silicon ! 10
ao
1.575 0.0006696284 4.271 Al
26.982 # octahedral aluminum ! 11
at
1.575 0.0009267943 3.302 Al
26.982 # tetrahedral aluminum ! 12
mgo
1.36 0.0004547008 5.264 Mg
24.305 # octahedral magnesium ! 13
mgh
1.05 0.0004547008 5.264 Mg
24.305 # hydroxide magnesium ! 14
cao
1.36
0.002532785 5.567 Ca
40.078 # octahedral calcium ! 15
cah
1.05
0.002532785 5.562 Ca
40.078 # hydroxide calcium ! 16
feo
1.575
0.004547008 4.906 Fe
55.845 # octahedral iron ! 17
lio
0.525
0.004547008 4.210 Li
S13
6.941 # octahedral lithium ! 18
Na
1
65.5126 2.350 Na
K
1
50.35557 3.334
Cs
1
50.35557 3.831 Cs
Ca
2
50.35557 2.872 Ca
Ba
2
23.66711 3.817 Ba
Cl
-1
50.40593 4.400 Cl
K
22.990 # aqueous sodium ion ! 19 39.098 # aqueous potassium ion ! 20 132.905 # aqueous cesium ion ! 21 40.078 # aqueous calcium ion ! 22 137.327 # aqueous barium ion ! 23 35.453 # aqueous chloride ion ! 24
! BOND: !! harmonic: K * (R - R0)**2 !!! K (float, Ken/angstrom^2): force constant !!! R0 (float, angstrom): equilibrium bond length !! fixed: R = R0 !!! R0 (float, angstrom): equilibrium bond length BOND
! 3
h*
o* harmonic
279038
1.000
! 1
ho
oh harmonic
279038
1.000
! 2
ho
ohs harmonic
279038
1.000
! 3
# ANGLE # h*
o*
h* harmonic 23048 109.47
# ao
oh
ho harmonic 15107 109.47
# ao
ohs ho harmonic 15107 109.47
# at
oh
# at
ohs ho harmonic 15107 109.47
# mgh oh
ho harmonic 15107 109.47
ho harmonic 15107 109.47
# mgh ohs ho harmonic 15107 109.47 # cah oh
ho harmonic 15107 109.47
S14
0.794
# cah ohs ho harmonic 15107 109.47 # feo oh
ho harmonic 15107 109.47
# feo ohs ho harmonic 15107 109.47 # lio oh
ho harmonic 15107 109.47
# lio ohs ho harmonic 15107 109.47 MOLECULE kaolinite
! SETTING: !! unitcell: unit cell dimensions (for crystalline materials only) !!! a (float, angstrom): unit cell edge length !!! b (float, angstrom): unit cell edge length !!! c (float, angstrom): unit cell edge length !!! alpha (float, degree): unit cell angle !!! beta (float, degree): unit cell angle !!! gamma (float, degree): unit cell angle !!! space_group (str, P1): space group SETTING
! 1
unitcell
5.154
8.942
7.401
91.69 104.61
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate ATOM name type x y z ! 34
S15
89.82 P1
Al1
ao
0.6650
4.3300
3.4040
! 1
Al2
ao
3.2140
2.8550
3.3960
! 2
Si1
st
4.9750
3.0050
0.6610
! 3
Si2
st
2.4620
1.4720
0.6710
! 4
O1
ob
-0.3210
3.0970
2.2630
! 5
O2
ob
0.0550
5.8590
2.2660
! 6
O3
ob
0.0140
4.4710
0.0000
! 7
O4
ob
1.0450
2.0680
0.1750
! 8
O5
ob
1.0710
6.8310
0.0020
! 9
O6
oh
-0.3200
8.5930
2.3290
! 10
O7
oh
3.8190
1.3530
4.3260
! 11
O8
oh
-0.9220
4.1030
4.3240
! 12
O9
oh
-0.9230
7.5290
4.3520
! 13
H1
ho
0.0860
0.2420
2.4870
! 14
H2
ho
-1.0150
1.4610
5.0180
! 15
H3
ho
-1.1230
4.1950
5.0680
! 16
H4
ho
-1.1110
6.9610
4.9970
! 17
Al3
ao
3.2560
8.8010
3.4040
! 18
Al4
ao
0.6510
7.3260
3.3960
! 19
Si3
st
2.4120
7.4760
0.6610
! 20
Si4
st
-0.1010
5.9430
0.6710
! 21
O10
ob
2.2700
7.5680
2.2630
! 22
O11
ob
2.6180
1.3880
2.2660
! 23
O12
ob
2.5770
0.0000
0.0000
! 24
O13
ob
3.6360
6.5390
0.1750
! 25
O14
ob
3.6340
2.3600
0.0020
! 26
O15
oh
2.2430
4.1220
2.3290
! 27
S16
O16
oh
1.2560
5.8240
4.3260
! 28
O17
oh
1.6690
8.5740
4.3240
! 29
O18
oh
1.6400
3.0580
4.3520
! 30
H5
ho
2.6770
4.7130
2.4870
! 31
H6
ho
1.5760
5.9320
5.0180
! 32
H7
ho
1.4690
8.6660
5.0680
! 33
H8
ho
1.4520
2.4900
4.9970
! 34
! BOND: BOND
! 8
H1
O6
! 1
0.750
H2
O7
! 2
0.770
H3
O8
! 3
0.776
H4
O9
! 4
0.880
H5
O15
! 5
0.750
H6
O16
! 6
0.770
H7
O17
! 7
0.776
H8
O18
! 8
0.880
# for use with CLAYFF MOLECULE Na
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate
S17
!! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate ATOM name type x y z ! 1 Na
Na
0.0000
0.0000
0.0000
! 1
# for use with CLAYFF MOLECULE Ca
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate !! z (float, angstrom): atom coordinate ATOM name type x y z ! 1 Ca
Ca
0.0000
0.0000
0.0000
! 1
# for use with CLAYFF MOLECULE Cl
! ATOM: !! name (str): atom name (required) !! type (str): atom type name !! charge (float, e): partial atomic charge !! x (float, angstrom): atom coordinate !! y (float, angstrom): atom coordinate
S18
!! z (float, angstrom): atom coordinate ATOM name type x y z ! 1 Cl
Cl
0.0000
0.0000
0.0000
! 1
Example LAMMPS Scripts System Setup # read lammps data file, write a restart file
# these variables must be set from the command line: variable
outname index init
variable
infile index init.data
# optional # initial temperature (do nothing if < 0) variable
beg_temp index -1
# these variables have generally conservative defaults: # kspace: 1 = pppm, 2 = ewald, 3 = ewald/n, other = off variable
kspace index 1
# restart != 0: write_restart at end of script variable
restart index 1
# ___end command line variable section___
# set up simulation # the following information is saved to restart files S19
units
real
atom_style
full
boundary
p p f
pair_style
lj/cut/coul/long 12
pair_modify
tail no
pair_modify mix arithmetic special_bonds lj 0.000000 0.000000 1.000000 special_bonds coul 0.000000 0.000000 1.000000 bond_style harmonic angle_style harmonic
read_data
${infile}
# end information that is saved to restart files
if "${kspace} == 1" then "kspace_style pppm 0.0001" if "${kspace} == 2" then "kspace_style ewald 0.0001" if "${kspace} == 3" then "kspace_style ewald/n 0.0001" kspace_modify slab 3.0
if "${beg_temp} >= 0" then & "velocity all create ${beg_temp} 24601 mom yes rot yes dist gaussian"
if "${restart} != 0" then & "write_restart
${outname}.restart.*"
S20
System Relaxation # NVE, NVT, NPT, or NPH ensemble
# these variables must be set variable outname index init.relax variable infile index init.restart.0 # ensemble/integrator: 1 = NVE, 2 = NVT, 3 = NPT, 4 = NPH, 5 = NVT/relax variable ensemble index 5 variable numsteps index 1000000
# damping terms are in multiples of ${timestep} # thermostat settings (if applicable) variable beg_temp index 1 variable end_temp index 320 variable tdamp index 10 # barostat settings (if applicable) variable beg_press index 200 variable
end_press index ${beg_press}
variable
pdamp index 1000
# ptensor = iso | aniso | tri | x | y | z | xy | xz | yz variable
ptensor index iso
# output options: # thermo_freq: output thermo snapshots at specified interval variable
thermo_freq index 1000
# dump_freq: dump trajectory snapshots at specified interval variable dump_freq index 1000 S21
# dump_detail: 1 = positions, 2 = +velocities, 3 = +forces variable
dump_detail index 1
# press_freq > 0: output pressure tensor data at specified interval variable
press_freq index 0
# other options: # kspace: 1 = pppm, 2 = ewald, 3 = ewald/n, other = off # (e.g set kspace=0 for non-periodic systems) variable
kspace index 1
# timestep: aka dt, e.g. 1 fs variable
timestep index 1
# processor allocation (i.e. "processors ${procx} ${procy} ${procz}") variable
procx index *
variable
procy index *
variable
procz index *
# restart_freq: save restart file every N steps (0 = only at end, = 0: reset current step to new_step variable
new_step index -1
# density > 0: rescale volume to result in specified density (avoid NPT, NPH) # note for ’real’ units: g/cm^3 * 0.6022 = g/mol/A^3 variable
density index 0
# ___end command line variable section___
# set up simulation
S22
processors ${procx} ${procy} ${procz} read_restart
${infile}
if "${kspace} == 1" then "kspace_style pppm 0.0001" if "${kspace} == 2" then "kspace_style ewald 0.0001" if "${kspace} == 3" then "kspace_style ewald/n 0.0001" neighbor
2.0 bin
neigh_modify
delay 5 check yes
timestep
${timestep}
kspace_modify slab 3.0
# change this to hold some atoms rigid #group MOBILE union all region MOBILE block INF INF INF INF 148 162 side out units box group MOBILE region MOBILE group FIXED subtract all MOBILE neigh_modify exclude group FIXED FIXED delete_bonds FIXED multi remove
# initial wall position should match previous final wall position variable zwall index "zlo 10" fix Zwall MOBILE wall/harmonic ${zwall} 50 0 1 units box extpress ${beg_press} 1e-8 variable FtoP equal 68568.411/lx/ly
# atm / (kcal/mol / Angstrom)
variable Pwall equal f_Zwall[1]*${FtoP}
# reset current simulation timestep if "${new_step} >= 0" then "reset_timestep ${new_step}"
S23
# pick desired ensemble/integrator variable _tdamp equal ${tdamp}*${timestep} variable _pdamp equal ${pdamp}*${timestep} if "${ensemble} == 1" then & "fix NVE MOBILE nve" & elif "${ensemble} == 2" & "fix TFIX MOBILE nvt temp ${beg_temp} ${end_temp} ${_tdamp}" & elif "${ensemble} == 3" &
"fix TFIX MOBILE npt temp ${beg_temp} ${end_temp} ${_tdamp} ${ptensor} ${beg_press} ${ elif "${ensemble} == 4" &
"fix TFIX MOBILE nph iso ${beg_press} ${end_press} ${_pdamp} nreset 1000 dilate partia elif "${ensemble} == 5" & "fix NVE MOBILE nve/limit 0.1" & "fix TFIX MOBILE langevin ${beg_temp} ${end_temp} ${_tdamp} 24601 tally yes" & else & "print ’Invalid ensemble/integrator!!!’" "exit"
# set up thermo data output if "${ensemble} > 1" then & "thermo_style custom step temp epair emol etotal press vol f_TFIX v_Pwall f_Zwall[2]" & else & "thermo_style custom step temp epair emol etotal press vol" thermo ${thermo_freq} thermo_modify flush yes
# set up dump/trajectory output -- potentially huge file!
S24
if "${dump_detail} == 3" then &
"dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu vx vy vz fx if "${dump_detail} == 2" then &
"dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu vx vy vz ma if "${dump_detail} == 1" then & "dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu mass q" dump_modify DUMP flush yes
# pressure tensor if "${press_freq} 0" then "restart ${restart_freq} ${outname}.restart" run
${numsteps}
if "${restart_freq} >= 0" then "write_restart ${outname}.restart.*"
NPT Simulation # NVE, NVT, NPT, or NPH ensemble
# naming conventions variable basename index t_330-p_20 variable suffix index 0 variable
outname index ${basename}.${suffix}
variable
infile index ${basename}.restart.*
# ensemble/integrator: 1 = NVE, 2 = NVT, 3 = NPT, 4 = NPH, 5 = NVT/relax variable ensemble index 2 variable numsteps index 7000000
# damping terms are in multiples of ${timestep} # thermostat settings (if applicable) variable beg_temp index 330 variable
end_temp index ${beg_temp}
variable
tdamp index 100
# barostat settings (if applicable) variable beg_press index 197.38465 variable
end_press index ${beg_press} S26
variable
pdamp index 1000
# ptensor = iso | aniso | tri | x | y | z | xy | xz | yz variable
ptensor index iso
# initial wall position should match previous final wall position variable zwall index "zlo 11.561567"
# output options: # thermo_freq: output thermo snapshots at specified interval variable
thermo_freq index 10000
# dump_freq: dump trajectory snapshots at specified interval variable
dump_freq index ${thermo_freq}
# dump_detail: 1 = positions, 2 = +velocities, 3 = +forces variable
dump_detail index 1
# press_freq > 0: output pressure tensor data at specified interval variable
press_freq index 0
# other options: # kspace: 1 = pppm, 2 = ewald, 3 = ewald/n, other = off # (e.g set kspace=0 for non-periodic systems) variable
kspace index 1
# timestep: aka dt, e.g. 1 fs variable
timestep index 1
# processor allocation (i.e. "processors ${procx} ${procy} ${procz}") variable
procx index *
variable
procy index *
variable
procz index *
# restart_freq: save restart file every N steps (0 = only at end, = 0: reset current step to new_step variable
new_step index -1
# density > 0: rescale volume to result in specified density (avoid NPT, NPH) # note for ’real’ units: g/cm^3 * 0.6022 = g/mol/A^3 variable
density index 0
# ___end command line variable section___
# set up simulation processors ${procx} ${procy} ${procz} read_restart
${infile}
if "${kspace} == 1" then "kspace_style pppm 0.0001" if "${kspace} == 2" then "kspace_style ewald 0.0001" if "${kspace} == 3" then "kspace_style ewald/n 0.0001" neighbor
2.0 bin
neigh_modify
delay 5 check yes
timestep
${timestep}
#change_box all z final -10 172 units box kspace_modify slab 3.0 diff ad
# change this to hold some atoms rigid #group MOBILE union all region MOBILE block INF INF INF INF 148 162 side out units box group MOBILE region MOBILE
S28
group FIXED subtract all MOBILE neigh_modify exclude group FIXED FIXED delete_bonds FIXED multi remove
fix Zwall MOBILE wall/harmonic ${zwall} 50 0 1 units box extpress ${beg_press} 1e-7 variable FtoP equal 68568.411/lx/ly
# atm / (kcal/mol / Angstrom)
variable Pwall equal f_Zwall[1]*${FtoP}
# reset current simulation timestep if "${new_step} >= 0" then "reset_timestep ${new_step}"
# pick desired ensemble/integrator variable _tdamp equal ${tdamp}*${timestep} variable _pdamp equal ${pdamp}*${timestep} if "${ensemble} == 1" then & "fix NVE MOBILE nve" & "compute TFIX_temp MOBILE temp" & elif "${ensemble} == 2" & "fix TFIX MOBILE nvt temp ${beg_temp} ${end_temp} ${_tdamp}" & elif "${ensemble} == 3" &
"fix TFIX MOBILE npt temp ${beg_temp} ${end_temp} ${_tdamp} ${ptensor} ${beg_press} ${ elif "${ensemble} == 4" &
"fix TFIX MOBILE nph iso ${beg_press} ${end_press} ${_pdamp} nreset 1000 dilate partia elif "${ensemble} == 5" & "fix NVE MOBILE nve/limit 0.1" & "fix TFIX MOBILE langevin ${beg_temp} ${end_temp} ${_tdamp} 24601 tally yes" & "compute TFIX_temp MOBILE temp" &
S29
else & "print ’Invalid ensemble/integrator!!!’" "exit"
# set up thermo data output if "${ensemble} > 1" then & "thermo_style custom step temp epair emol etotal press vol f_TFIX v_Pwall f_Zwall[2]" & "thermo_modify temp TFIX_temp" & else & "thermo_style custom step temp epair emol etotal press vol" thermo ${thermo_freq} thermo_modify flush yes
# set up dump/trajectory output -- potentially huge file! if "${dump_detail} == 3" then &
"dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu vx vy vz fx if "${dump_detail} == 2" then &
"dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu vx vy vz ma if "${dump_detail} == 1" then & "dump DUMP all custom ${dump_freq} ${outname}.lammpstrj id type mol xu yu zu mass q" dump_modify DUMP flush yes
# pressure tensor if "${press_freq} 0" then "restart ${restart_freq} ${basename}.restart.*" run
${numsteps}
if "${restart_freq} >= 0" then "write_restart ${basename}.restart.*"
Acknowledgement This material is based upon work supported as part of the Center for Frontiers of Subsurface Energy Security, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0001114. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DEAC04-94AL85000.
S31
References (1) Neder, R. B.; Burghammer, M.; Grasl, T.; Schulz, H.; Bram, A.; Fiedler, S. Refinement of the Kaolinite Structure from Single-Crystal Synchrotron Data. Clays Clay Miner. 1999, 47, 487–494. (2) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. Molecular Models of Hydroxide, Oxyhydroxide, and Clay Phases and the Development of a General Force Field. J. Phys. Chem. B 2004, 108, 1255–1266. (3) Cygan, R. T.; Greathouse, J. A.; Heinz, H.; Kalinichev, A. G. Molecular Models and Simulations of Layered Materials. J. Mater. Chem. 2009, 19, 2470–2481. (4) Thyveetil, M.-A.; Coveney, P. V.; Greenwell, H. C.; Suter, J. L. Computer Simulation Study of the Structural Stability and Materials Properties of DNA-Intercalated Layered Double Hydroxides. J. Am. Chem. Soc. 2008, 130, 4742–4756. (5) Wang, J. W.; Kalinichev, A. G.; Kirkpatrick, R. J.; Cygan, R. T. Structure, Energetics, and Dynamics of Water Adsorbed on the Muscovite (001) Surface: A Molecular Dynamics Simulation. J. Phys. Chem. B 2005, 109, 15893–15905. (6) Greathouse, J. A.; Cygan, R. T. Molecular Dynamics Simulation of Uranyl(VI) Adsorption Equilibria onto an External Montmorillonite Surface. Phys. Chem. Chem. Phys. 2005, 7, 3580–3586. (7) Berendsen, H.; Postma, J.; Van Gunsteren, W.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermolecular forces 1981, 11, 331–342. (8) Teleman, O.; Jonsson, B.; Engstrom, S. A Molecular-Dynamics Simulation of a Water Model with Intramolecular Degrees of Freedom. Mol. Phys. 1987, 60, 193–203. (9) Raabe, G.; Sadus, R. J. Molecular Dynamics Simulation of the Dielectric Constant of Water: the Effect of Bond Flexibility. J. Chem. Phys. 2011, 134, 6. S32
(10) Raabe, G.; Sadus, R. J. Molecular Dynamics Simulation of the Effect of Bond Flexibility on the Transport Properties of Water. J. Chem. Phys. 2012, 137, 8. (11) Wallqvist, A.; Teleman, O. Properties of Flexible Water Models. Mol. Phys. 1991, 74, 515–533. (12) Smith, D. E.; Haymet, A. D. J. Simulations of Aqueous-Solutions - the Role of Flexibility and the Treatment of Long-Range Forces. Fluid Phase Equilib. 1993, 88, 79–87. (13) Koneshan, S.; Rasaiah, J. C.; Lynden-Bell, R. M.; Lee, S. H. Solvent Structure, Dynamics, and Ion Mobility in Aqueous Solutions at 25 Degrees C. J. Phys. Chem. B 1998, 102, 4193–4204. (14) Smith, D. E.; Dang, L. X. Computer-Simulations of Nacl Association in Polarizable Water. J. Chem. Phys. 1994, 100, 3757–3766. (15) Cygan, R. T.; Romanov, V. N.; Myshakin, E. M. Molecular Simulation of Carbon Dioxide Capture by Montmorillonite Using an Accurate and Flexible Force Field. J. Phys. Chem. C 2012, 116, 13079–13091. (16) Zhu, A.; Zhang, X.; Liu, Q.; Zhang, Q. A Fully Flexible Potential Model for Carbon Dioxide. Chin. J. Chem. Eng. 2009, 17, 268–272. (17) Harris, J. G.; Yung, K. H. Carbon Dioxides Liquid-Vapor Coexistence Curve and Critical Properties As Predicted by a Simple Molecular-Model. J. Phys. Chem. 1995, 99, 12021–12024. (18) Plimpton, S. J. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. (19) Plimpton, S. J. LAMMPS Molecular Dynamics Simulator. http://lammps.sandia. gov.
S33
(20) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157–2164.
S34
