Jul 16, 2015 - convergence criteria (energy change per atom ...

0 downloads
0 Views
4MB Size

This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.

Article pubs.acs.org/JPCC

Systematic Optimization of a Force Field for Classical Simulations of TiO2−Water Interfaces Erik G. Brandt and Alexander P. Lyubartsev* Department of Materials and Environmental Chemistry, Stockholm University, SE 106 91, Stockholm, Sweden S Supporting Information *

ABSTRACT: Atomistic force ﬁeld parameters were developed for the TiO2−water interface by systematic optimization with respect to experimentally determined crystal structures (lattice parameters) and surface thermodynamics (water adsorption enthalpy). Optimized force ﬁeld parameters were determined for the two cases where TiO2 was modeled with or without covalent bonding. The nonbonded TiO2 model can be used to simulate diﬀerent TiO2 phases, while the bonded TiO2 model is particularly useful for simulations of nanosized TiO2 and biomatter, including protein− surface and nanoparticle−biomembrane simulations. The procedure is easily generalized to parametrize interactions between other inorganic surfaces and biomolecules.

■

generalized to develop force ﬁeld parameters for other inorganic surfaces and biomolecules. In an accompanying paper,13 we use the optimized force ﬁeld parameters to address the other long-standing problem in nanobio simulationsaccurate sampling of biomolecules near inorganic surfacesand to compute binding free energies of amino acid side chain analogues and a peptide to the TiO2 surface. TiO2 is one of the most common materials forming interfaces with biomatter, and is found in consumer products such as paints, sunscreens, and food coloring.14 Bulk TiO2 is believed to be inert with respect to biomaterials, but its role in nanotoxicity is still unclear. There has been considerable interest in developing models for TiO2 under ambient conditions, because of its abundant presence in biological environments. The most widely used classical model is due to Matsui and Akaogi15 (MA), which was parametrized to reproduce the bulk structure of TiO2 (not any surface properties). Later eﬀorts extended the MA model with water interactions16−18 but were based on pairwise models that are diﬃcult to unite with the standard force ﬁelds for biomolecular simulations already available in the literature (e.g., AMBER,19 CHARMM,20 OPLS21), and/or requires new parametrizations for every additional atom type that is to be simulated together with TiO2. Other simulations of the TiO2/water interface combined force ﬁeld parameters from solid (TiO2) and liquid (water) state, in which case solid bulk atoms had to be ﬁxed to keep the bulk structure of TiO2 intact.22 More complex TiO2 force ﬁelds oﬀer the possibility to model surface reactions23 and adsorption24 but at the cost of

INTRODUCTION Biomolecules interacting with inorganic objects are fundamental in nanobiotechnological applications such as surfaceattached biomolecules for target delivery of drug molecules,1 protein adsorption to medical implants,2 protein−nanoparticle interactions3 and, not least, to understand the molecular mechanisms of nanotoxicity.4 The contact points between inorganic surfaces and biomatter (the nanobio interface5) can be tracked with experimental techniques like dynamic light scattering (DLS),6−8 chromatography9,10 and/or spectroscopy,11,12 but only by indirect measurements. On the other hand, computer simulations can in principle be used to directly calculate interactions at the nanobio interface, covering system sizes ∼1000 Å and time scales ∼1000 ns. These windows are highly relevant when modeling the nanobio interface, but outside the realm of quantum mechanics. This calls for “classical” models with atomistic or semiatomistic (atoms being grouped into eﬀective interaction centers) representations in computer simulations of nanobio interactions, parametrized to reproduce properties relevant to the nanobio interface. Substantial eﬀort has been invested during the last decades to develop classical molecular models for biomolecules (proteins, lipids, nucleic acids, carbohydrates, etc.), but there has been no comparable endeavor to include inorganic materials into the models. The state of classical models aimed to describe the nanobio interface remains underdeveloped, largely because of the diﬃculties involved with representing the electronic structure of the inorganic material by interacting point particles. The present work is an attempt to close this modeling gap by developing atomistic force ﬁeld parameters for the TiO2−water interface. We use an automated algorithm that is extended to accept arbitrary experimental data as targets in the parameter ﬁtting. The method can easily be © 2015 American Chemical Society

Received: March 19, 2015 Revised: June 27, 2015 Published: July 16, 2015 18110

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

i.e., parametrized at speciﬁc temperatures and pressures, and not necessarily accurate at other state points. The energy expression in eq 1 has been chosen to ensure transferability with the widespread AMBER force ﬁeld,19 which includes interactions parameters for lipids, proteins, nucleic acids, and other common organic molecules. New compounds can also be added to AMBER with relatively small eﬀort. We note that the energy functions for the other common force ﬁelds in the simulation community (CHARMM,20 OPLS21) are very similar. Finally, biological relevance and experimental basis were met by parametrizing the TiO2 force ﬁeld to the adsorption enthalpy of water at room temperature and normal pressures. DFT Geometry Optimizations of Water Molecules on TiO2 Surfaces. The Materials Studio27 (MS) program was used to prepare and run the periodic density functional theory (DFT) calculations. The rutile TiO2 unit cell belongs to space group P42/mnm, with lattice parameters a = b = 4.594 Å and c = 2.959 Å (the unit cell angles are all 90°). The unit cell was replicated and cleaved either along the (100) plane or along the (110) plane. Water molecules were placed in random (unfavorable) orientations close to the surface. Periodicity was enforced in all dimensions, and the simulation cell was elongated to 15 Å in the direction normal to the surfaces. The sizes of the ﬁnal TiO2 slabs were 9.1874 × 5.9174 Å2 (100 surface) and 5.9174 × 6.4965 Å2 (110 surface). The CASTEP28 module as implemented in MS was used for geometry optimizations with DFT using the Perdew−Wang 91 exchange-correlation functional29 and ultrasoft pseudopotentials in reciprocal space representation. There were in total 200 electrons in the (100) system and 176 electrons in the (110) system. The cell size and the positions of the surface atoms were held ﬁxed during the optimization to avoid excessive rearrangements of surface atoms. The optimal geometries were determined self-consistently by the BFGS solver in CASTEP, which took about 100 iterations with tight convergence criteria (energy change per atom <10−5 eV, the maximum force was <3 × 10−2 eV/ Å and the maximum displacement was <10−3 Å). TiO2 Bulk Models. TiO2 has been studied with computational methods including ab initio calculations16 and classical molecular dynamics (MD) simulations.24 The interaction parameters employed in MD studies have usually been determined from quantum mechanical (QM) or semiempirical calculations, and the energetic interactions in TiO2 have often been assumed to be completely noncovalent.15 However, the short Ti−O distance (1.96 Å) in the crystal structure, and the experimentally determined charge distribution of rutile30 demonstrates that there is a substantial electron overlap between the Ti and O atoms. The covalent contribution to the Ti−O bond is also corroborated by calculations based on semiempirical basicity/acidity measurements.31,32 Neglecting the covalent contribution leads to a model for the metal oxide with an iono-covalent imbalance. In a bulk system, this can be compensated by tuning the nonbonded parameters accordingly, but when an interface is present in the system, such tuning becomes impossible. We consider two models in this paper. Our primary concern is to parametrize a novel bonded TiO2 model with a systematic optimization procedure. In addition, we use the same method and show that the widespread Matsui−Akaogi (MA) potential can be improved for interface simulations by systematic

nonstandard functional forms or multibody potentials that may have limited transferability. The purpose of this paper is therefore to determine interaction parameters for classical MD simulations of TiO2 surfaces in contact with aqueous media (including biomatter) based on the following criteria: (1) compatibility with standard force ﬁelds for biomolecular simulations; (2) transferability of TiO2 interactions to a wide variety of biomolecules: proteins, lipids, and nucleic acids; nonbonded interactions with new molecules should be deduced by combination rules; (3) The force ﬁeld should reproduce interactions at ambient conditions, i.e., room temperature and atmospheric pressure, which are biologically relevant; and (4) The force ﬁeld parameters should be derived in a thermodynamically consistent manner, and based to the furthest extent on known experimental data. The modeling is underpinned by quantum mechanical geometry optimizations of water on TiO2 surfaces using density functional theory. A bonded TiO2−surface model is developed based on the optimized geometries. The model parametrization is carried out by a systematic approach,25,26 which we modify to target a combination of structural (lattice parameters) and thermodynamical (adsorption enthalpies) TiO2 data. The next section, “ Methods and Materials”, describes the basic modeling approach, the involved interaction parameters, and the method used to perform systematic optimization of the parameters. The “Results and Discussion” section presents the outcome of the optimizations, discusses their optimal values, and investigate the structure and dynamics of the TiO2− water interface for three diﬀerent TiO2 models. “Conclusions” summarizes the most important conclusions and outlines the next steps for further studies.

■

METHODS AND MATERIALS To fulﬁll the compatibility criteria listed in the introduction, we start from the standard molecular mechanical force ﬁeld expression:

∑

U=

kb , ij(rij − r0, ij)2

ij bonded

+

∑

∑

kθ , ijk(θijk − θ0, ijk)2 +

ijk bonded

ijkl bonded

1 Vφ , ijkl 2

∑

× [1 + cos(nφijkl − φ0, ijkl )] +

4ϵij

ij nonbonded(excl)

⎡⎛ σij ⎞12 ⎛ σij ⎞6 ⎤ × ⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥ + ⎢⎣⎝ rij ⎠ ⎝ rij ⎠ ⎥⎦

∑ ij nonbonded(excl)

1 qiqj 4π ϵ0 rij (1)

The ﬁrst three terms are bonded, angle, and dihedral interactions, while the last two terms are nonbonded: Lennard-Jones (12−6) and Coulomb interactions. The force ﬁeld consists of the sets of force constants (kb, kθ, and Vφ), equilibrium distances (r0, θ0 and φ0), nonbonded parameters (ϵ and σ) and partial charges (q) in eq 1. The major limitations to the scope of simulations based on eq 1 are as follows: (1) No modeling of chemical reactions involving bond formation/breaking because of predetermined covalent coordination of the atoms. The surface chemistry in terms of hydroxylation and protonation is ﬁxed prior to the simulation. (2) Electrostatic interactions are approximated with ﬁxed partial charges, ignoring charge transfer and including polarization only in an implicit fashion. (3) The force ﬁeld parameters in eq 1 are state point-dependent, 18111

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C optimization of its original parametrization. The original MA model is nonbonded and based on the Buckingham potential UBH(rij) = Aij exp( −rij/bij) − Cij/rij6

(2)

but we will use the Lennard-Jones potential form (cf. eq 1) ⎡⎛ σij ⎞12 ⎛ σij ⎞6 ⎤ ULJ(rij) = 4ϵij⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥ ⎢⎣⎝ rij ⎠ ⎝ rij ⎠ ⎥⎦

(3)

with ﬁtted parameters for Ti and O and combination rules to determine the cross-terms. The ﬁt is excellent and we consider the two variants of the MA model to be equivalent. The LJ parameters for the MA model are listed in Table 1. The partial Table 1. Buckingham and Lennard-Jones Parameters for the Matsui−Akaogi model15 atom

q (e)

A (kJ mol−1)

b (Å)

C (kJ mol−1)

σ (Å)

ϵ (kJ mol−1)

Ti O

2.196 −1.098

3002664 1136872

0.154 0.234

506 2916

1.958 2.875

2.5422 1.3897

Figure 1. Ti is 6-coordinated and O is 3-coordinated in the bulk. The coordination number for Ti is 4 or 5 on surfaces, while O is 2-coordinated.

charges in the MA model are +2.196e for Ti and −1.098e for O, originally determined from ﬁtting phonon dispersion curves from inelastic neutron scattering to a rigid ion model.33 The rest of the parameters were determined by ﬁtting simulation data to the experimentally determined lattice parameters and elastic constants (by hand). We note that, by neglecting the covalent contribution to the Ti−O bond artiﬁcially, small atom radii (σ’s) are needed to obtain the predicted TiO2 structure in the MA model (Table 1). The MA model can not be expected to perform well for interface simulations, since the original parametrization was done for bulk structures. In our bonded model, we adopt a ﬁxed (harmonic) Ti−O bond. No angular or dihedral restraints were added, so all kθ = Vφ = 0 in eq 1. To account for surface chemistry, we deﬁne atom types based on coordination and the following constraints on the atoms and the charges in the system. 1 The total ratio of O to Ti atoms in the system is exactly 2:1. 2 Ti atoms are 4-to-6-coordinated, and O atoms are 2-to-3coordinated, and there are two 2-coordinated O atoms per 4-coordinated Ti atom and one 2-coordinated O atom per 5-coordinated Ti atom (see Figure 1). 3 The charges in the bulk are related by q3 = −q6 /2

This gives us 5 atom types, 2 for O (OT2 and OT3) and 3 for Ti (Ti4, Ti5 and Ti6). We will perform the parametrization for the (100) surface where there are no 4-coordinated Ti and can ignore this atom type for the remainder of the paper. The parameters to be determined for the bonded model are then LJ parameters and charges (3 parameters per atom type), and bond parameters for the Ti−O bond (2 parameters, the same bond is used between all Ti and O atom types). The nonbonded interaction terms between Ti and O are determined by the standard Lorentz−Berthelot combination rules: σij = (σi + σj)/2 and ϵij = ϵiϵj . No Lennard-Jones interactions between atoms connected with less than three bonds are included in the energy calculation. Assigning physically motivated parameters to eq 1 yields a realistic covalent/ionic proportion that accurately reproduces surface interactions. We adopt the parametrization procedure outlined by Heinz et al.34 This strategy emphasizes the assignment of accurate partial atomic charges to the force ﬁeld (qi in eq 1), because Coulomb interactions outweigh Lennard-Jones interactions with roughly 2 orders-of-magnitude in inorganic polar materials. To a high degree, surface properties are determined by electrostatic interactions, while bonded parameters play the lesser role of enforcing the atomic structure within the bulk material. The ﬁnal force ﬁeld parameters are summarized in Table S1 in the Supporting Information for the reader in a hurry. TiO2 surfaces. Our aim is to parametrize TiO2/water interactions under ambient conditions, i.e., neutral pH and salt concentrations up to 150 mM. Titanium surfaces in aqueous media are expected to be oxygen-covered under such conditions,35−37 e.g., as implants38 or nanoparticle surfaces. Therefore, we consider an oxygen-terminated (100) TiO2 surface consisting of a layer of 2-coordinated oxygens on top of 5-coordinated titaniums (Ti−O−Ti bridges, similar to the setup in ref 22. The surface chemistry of titania is largely determined by the hydroxylation/-protonation of surface sites (surface pH) and the degree of surface ionization. Biological salt concentration (150 mM) corresponds to 1 salt molecule per 10 nm3, which

(4)

so that the basic TiO2 unit is neutral. 4 The system is neutral, 6

∑ qiNi = 0 i=2

(5)

i.e., the total charge in the system is zero. We can take q6 and q5 as the two charges to be varied, and combine the previous conditions to derive the following relations for the other charges: q2 =q6 /2 − q5 q3 =−q6 /2 q4 =2q5 − q6

(6) 18112

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

by ϵ. The exact values of these parameters depends on many factors, in particular the local geometry and coordination chemistry. Tabulated σ-values have been standard for the lighter elements.48,49 Alvarez50 recently reported revised vdW radii for the 93 naturally occurring elements in the periodic table based on a statistical analysis of the Cambridge Structural Database. We used Alvarez’s value, σTi = 4.92 Å (and σO = 3.0 Å, in excellent agreement with all atomistic protein force ﬁelds) as the starting value for our force ﬁeld. ϵ represents the atom polarizability and depends more sensitively than σ on local geometry. It is expected to be lower within the TiO2 covalent framework than, e.g., pure fcc elemental metals.51 A very rough initial estimate of ϵTi can be obtained from the London dispersion equation,42 Edisp = 4ϵσ6 = (3/4)α02I, where α0 is the electronic polarizability per volume (in units of 4πϵ0), and I is the ﬁrst ionization potential of the atom. The experimental52 values α0,Ti = 14.6 Å3 and ITi = 658.2 kJ/mol yield ϵTi = 1.85 kJ/mol using the vdW radius from ref 50. We note that this is only an order-ofmagnitude estimation, because the same calculation for oxygen (α0,O = 0.8 Å3 and IO = 1313 kJ/mol) yields ϵO = 0.216 kJ/mol, while ϵ is typically 3−4 times larger in organic compounds (e.g., ϵO = 0.636 kJ/mol in TIP3P water). Bond Parameters. The initial value for the equilibrium bond distance was chosen according to the crystal structure of rutile TiO2.53 The unit cell geometry is described by the rutile Ti−O distance, rTi−O = 1.96 Å(and the two angles ∠Ti−O−Ti = 90° and ∠O−Ti−O = 99°; see, e.g., Figure 2 in ref 54). The force constant for the Ti−O bond was obtained from the TiO2 IR spectrum,55 which is dominated by a peak at ν ∼ 1000 cm−1, corresponding to a frequency f = cν = 30 THz. This relates to the bond force constant in eq 1 by kb = μω2/2, where μ = mTimO/(mTi + mO) is the reduced mass and ω = 2 π f is the angular frequency. The IR frequency together with mTi = 48u and mTi = 16u (where u = 1.660 × 10−27 kg is the uniﬁed atomic mass unit) yields kb = 2131.00 kJ/mol/Å2. Systematic Optimization Procedure. The ﬁnal force ﬁeld parameters were determined by systematic optimization of the initial values just described. The simpler MA model was optimized by varying the parameters simultaneously in a single sweep, but a two-step procedure was employed for the bonded model: 1. Interaction parameters for bulk TiO2 were obtained by systematic optimization of simulated (rutile) crystal structures to experimentally determined lattice parameters. 2. The conﬁgurations obtained from quantum-mechanical geometry optimizations of water molecules on TiO2 slabs were used as starting structures in systematic optimization of TiO2−H2O interaction parameters. The aim of the optimization algorithm is to ﬁnd a set of M parameters λopt that best reproduce N reference measurements Ri, i = 1,..., N. This is achieved by constructing an objective function

translates to at most 3 molecules for the small water volumes employed here, and will be neglected for the remainder of this paper. The deﬁnition of pH in simulations is not straightforward due to the lack of free proton concentration. Ionizable groups are usually predetermined in protonated/deprotonated forms depending on the relationship between the chosen pH and pK of the ionizable group. Experimental validation of TiO2/water surface chemistry is provided by the ζ-potential, which measures the electrical potential diﬀerence between the solvent and the stationary layer of ﬂuid attached to the TiO2 surface (the slipping plane). The pH-value where the ζ-potential is zero is called the point of zero charge (PZC) and represents protonation of half of the surface sites. The TiO2 ζ-potential decreases with increasing pH and levels out around −50 mV at neutral pH 7 (e.g., Figure 7 in ref 39) after which there is no further surface protonation. The PZC is around pH 5, as conﬁrmed by previous experiments40 and theoretical calculations.41 A simple estimate of the fraction α of disassociated (deprotonated) surface sites at neutral pH is then42 α = 10−5/10−5 + [H+]) = 0.99, i.e., almost all surface sites are deprotonated. Theoretical modeling of hydroxylation/protonation of the (100) surface also corroborates this picture (see Figure 2B in ref 22). Hence, we modeled our surface without OH-groups to mimic neutral pH based on the following motivations. (1) Experimentally determined ζ-potential levels out at neutral pH 7 indicating that only a small number of surface oxygens are protonated. (2) Given the PZC for TiO2 around pH 5, 1% protonation at pH 7 is expected. (3) 2-coordinated oxygens with ﬁlled electron shells in Ti−O−Ti bridges on the surface are highly unlikely to be protonated, as demonstrated by ab initio molecular dynamics (AIMD) simulations, which calculated pK of such protonation to be −1.43 (4) We aim to mimic oxygen-rich conditions as to be expected in biological environments. Atomic Parameters. Atomic Charges. Surface properties of polar materials are dominated by electrostatic interactions and are modeled by the partial charges, qi, in eq 1. The accurate assignment of these parameters are therefore of high importance for the overall performance of the force ﬁeld with regards to surface properties. Within the harmonic approximation used in eq 1, the electrostatic interactions in bulk TiO2 are completely determined by the partial charge q6 (assuming charge neutrality q3 = −q6/2). Previous modeling eﬀorts based on quantum mechanical calculations have assigned a wide range of values to q6 (= qTi), from 0.68e44 to 2.2e.16 This highlights the diﬃculty for ab initio calculations to determine electron densities with accuracy that allows the extraction of robust partial charges for inorganic crystals, in particular for elements containing d- and f-electrons.45 The charge density of TiO2 rutile has been carefully determined by Jiang et al.30 using a combination of electron and X-ray diﬀraction. The electron deformation densities were ﬁtted to spherical atomic basins46 using multipole reﬁnement,47 where the optimal reﬁnement parameters furnish valence populations (i.e., partial charges). The multipole analysis yielded q6 = 1.38 ± 0.08 e for Ti and q3 = −0.69 ± 0.04 e30 for O. The uncertainties are due to a small part of the charge distribution, which can not be unequivocally distributed within a spherical model. We note that this is a substantial reduction from the formal Ti charge (+4e), and suggests a substantial covalent contribution to the Ti−O bond. This observation is conﬁrmed by semiempirical calculations based on oxide acidity/basicity scales.31,32 Lennard-Jones Parameters. The atomic size is represented by σ, while the strength of the dispersion energy is determined

⎛ R i − Si(λ) ⎞2 χ (λ) = ∑ ⎜ ⎟ + σi ⎠ i=1 ⎝ N

2

⎛ λ − λ ⎞2 ∑ ⎜⎜ m m,0 ⎟⎟ α λm ⎠ m=1 ⎝ M

(7)

and then solving the nonlinear minimization problem λopt = min χ 2 (λ) λ

(8)

by iteration. Si(λ) is the calculated value of the reference data from a simulation with the current (suboptimal) value of the parameter set 18113

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

Table 2. Target Data and Final Reproduction of Target Data after Systematic Optimization for the MA Model and the Bonded TiO2 Model target

units

experimental

scaling factor

Matsui−Akaogi

nonbonded

bonded

a b c ΔHads

Å Å Å kJ mol−1

4.592 4.592 2.958 95.6

0.04592 0.04592 0.02958 5

+2.1% +2.1% +6.6% +72%

−1.2% −1.2% +3.8% +15%

+1.2% +1.2% −3.0% −3.1%

λ. The second term in eq 7 is needed to avoid overﬁtting and puts a harmonic penalty for parameters that stray too far from their initial values. αλm is the prior of parameter λm and deﬁnes an interval around its initial value λm,0 in which the parameter is expected to vary. The priors are provided before the optimization and are selected by hand from chemical and physical knowledge of the system. We will refer to the reference data Ri as targets, which will be described in detail shortly. Each term in the objective function eq 7) is scaled with the factor σi (Table 2)to give the targets similar magnitudes and prevent ill-conditioned matrix inversion during the optimization. The scaling factor is determined from the uncertainty in the measurement of the target data. The minimization problem stated in eq 8 was solved with ForceBalance,25,26 which is open-source software that provides a framework for force ﬁeld optimization. ForceBalance has been designed to accept user-deﬁned targets and combine reference data obtained by experimental measurements and ab initio calculations. The χ2-minimization can be performed by any suitable algorithm, e.g.,ﬂavors of quasi-Newton or steepest descent. The major technical diﬃculty in practice is evaluating the numerical derivatives ∇λSi of the simulation data with respect to the force ﬁeld parameter λn, since Si is subject to thermal noise. When Si is the ensemble average of a quantity f(r, V), which not explicitly depend on λ, its analytic derivative with respect to λ (where we have suppressed the subscript n) can be written56 d ⟨f ⟩ Si, λ = dλ d ⎧ 1 ⎨ = dλ ⎩ Q (λ) ⎡ dU = −β ⎢ f ⎣ dλ

We extended ForceBalance with an in-house written module to accept a general target type for thermodynamic ﬁtting. It is based on analytical derivatives such as eq 10 for force ﬁeld optimization to thermodynamic quantities. Further details can be found in the Supporting Information. Target Data. Both TiO2 models (nonbonded and bonded) were optimized to the experimentally determined lattice parameters of rutile57 and to the adsorption enthalpy of water on rutile at 300 K58 (see Table 2). The rutile crystal unit is rectangular with a = b = 4.592 Å and c = 2.958 Å. The lattice parameters were extracted from the box vectors in molecular dynamics simulations at constant pressure of rutile (5 × 5 × 8)supercells (see “ Simulation Protocol for Optimization section”). The TiO2−water enthalpy measures the enthalpy diﬀerence for water molecules at ﬁxed concentration adsorbed to the TiO2 surface compared to the gas phase. The dilute limit represents the enthalpy diﬀerence for a single water molecule on the TiO2 surface as measured in microcalometric experiments.58 The water adsorption enthalpy is extracted from three independent simulations: (1) An interface simulation of the TiO2 slab and the water molecule, (2) a surface simulation of the isolated TiO2 slab, and (3) a gas simulation of a single water molecule (which could also be estimated from an ideal gas calculation). The adsorption enthalpy is obtained as Δhads = Hadsorbed /Nwater − hgas = (Hinterface − Hsurface)/Nwater − hwater

∫

⎫ dr dVf (r, V ) exp[− βU (r, V ; λ) − βpV ]⎬ ⎭ dU ⎤ − ⟨f ⟩ ⎥ dλ ⎦ (9)

where Nwater = 1 is the number of water molecules in the simulation (and hwater is the enthalpy of a gas molecule). It is assumed that the number of TiO2 units in the solid slab are equal in the interface and surface simulations. Analytic derivatives for some thermodynamic properties (including density and heat of vaporization) were given in ref 26, but not for lattice parameters or adsorption enthalpy. Let f(r,L) = L be any box vector in eq 9, in which case

where the brackets in the last line denotes averages in an ensemble with ﬁxed λ-value. Q(λ) is the canonical partition function, and U(r, V; λ) is the total potential energy function (eq 1), and the volume-dependences of U and f are included to account for simulations with periodic boundary conditions. β = 1/kBT (the inverse of the product of the Boltzmann constant and the absolute temperature) and the applied external pressure p are ﬁxed control parameters that deﬁne the ensemble. In other words, the derivative Si,λ is the diﬀerence between two ensemble averages and can be evaluated without imprecise numerical diﬀerences. When f is a linear function of U, i.e., f(r,V;λ) = aU(r,V;λ) + bg(r,V), where a and b are constants and g(r,V) is a function that does not explicitly depend on λ (but f does), a similar formula applies, according to fλ = a

dU dλ

⎡⎛ dU − β ⎢⎜a U ⎝ dλ ⎣

⎛ dU − ⎜a⟨U ⟩ ⎝ dλ

+ b⟨g ⟩

+b g d U ⎞⎤ ⎟⎥ d λ ⎠⎦

(11)

⎡ dU S1, λ = Lλ = −β ⎢ L ⎣ dλ

− ⟨L⟩

dU ⎤ ⎥ dλ ⎦

(12)

for the lattice parameters. Since the adsorption enthalpy is a combination of three contributions and a linear function of the potential energy, H = U + pV, the expression for the derivative of h is more complicated. Combining eq 10 and eq 11 leads to dUi dU dUi ⎤ β⎡ − ⎢ Ui i − ⟨Ui⟩ ⎥ N⎣ dλ dλ dλ ⎦ dU dUs ⎤ β⎡ 1 dUs − + ⎢ Us s − ⟨Us⟩ ⎥ N dλ N⎣ dλ dλ ⎦ dUs dUs ⎤ βp ⎡ βp ⎡ dUi dUi ⎤ − − ⟨Vi ⟩ − ⟨Vs⟩ ⎢ Vs ⎥ ⎥+ ⎢ Vi ⎦ ⎣ ⎣ N N dλ dλ dλ dλ ⎦ ⎡ dUg dUg dUg ⎤ ⎥ − + β ⎢ Ug − ⟨Ug⟩ ⎢⎣ dλ dλ dλ ⎥⎦ (13)

S2, λ = Δhads, λ =

dU ⎞ ⎟ dλ ⎠ (10) 18114

1 N

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

what decisions we therefore made in constructing the models. We then show that our method can be used to improve the MA model substantially for interface simulations without any additional parameters. Finally, we report on the full optimization of our bonded TiO2 model. Optimal Geometries of Water Molecule on Rutile Surfaces. We performed the geometry optimizations for a H2O molecule on two diﬀerent TiO2 surfaces with the aim to ﬁnd how the adsorption depends on whether oxygen or titanium atoms are exposed at the surface. We chose the (100) and the (110) surfaces of TiO2 to represent the two cases. The (100) surface have 2-coordinated oxygen atoms and 5-coordinated titanium atoms at the surface. The (110) surface is a more complex combination of dangling titanium and oxygen atoms with diﬀering coordination numbers. The optimized geometries of the water molecules on the surfaces, as obtained from DFT calculations, are shown in Figure 3a

for the adsorption enthalpy, where we have used the shorthand notations N = Nwater, and i = interface, s = surface and g = gas for the diﬀerent simulations. Simulation Protocol for Optimization. Each simulation was performed on a 8-core machine using the GROMACS 4.6.5 package59,60 as engine. The VESTA program61 was used to create TiO2 crystal and slab structures from replication of the rutile unit cell.55 A 5 × 5 × 8-supercell (with periodic atoms at the boundaries removed) was used to extract accurate lattice parameters from bulk simulations (Figure 2(a)). The adsorption

Figure 2. Systems used to calculate target data (lattice parameters and TiO2−water adsorption enthalpy) during the optimization of the TiO2 force ﬁelds. (a) Lattice parameters were extracted from bulk simulations of a (5 × 5 × 8)-supercell of rutile TiO2 with periodic boundaries. (b) TiO2−water adsorption enthalpy was extracted from interface simulations of a (5 × 8) TiO2 slab.

enthalpy was calculated from a thick (100) slab (Figure 2b), with and without a water molecule, according to eq 11. The (100) surface is terminated in both a1-directions with oxide layers. The slab is symmetric, consists of stacked Ti and O layers along axis a1 (normal to the surface), and is periodic in the lateral directions. Covalent bonds were assigned to Ti and O atoms within 2 Å. Bulk Ti and O atoms were then 6- and 3-coordinated, respectively, and the corresponding surface atoms were 5- and 2-coordinated. The systems were equilibrated for 20 ps (bulk simulations) or 100 ps (interface simulations) at the current values of the varied force ﬁeld parameters, before target data was calculated from production simulations lasting 50 ps (bulk) or 100 ps (interface). All simulations were performed at constant temperature (300 K) with the thermostat of Bussi et al.62 and constant anisotropic pressure tensor (NPT ensemble) using Berendsen’s weak scaling algorithm.63 All the elements of the reference pressure tensor were set to 1 bar, while the diagonal elements of the compressibility tensor were set to 5 × 10−7 bar−1 (corresponding to the inverse of the TiO2 bulk modulus KTiO2 = 200 GPa64). The oﬀ-diagonal elements were set to zero to ensure a rectangular simulation box (tests with an anisotropic compressibility tensor showed no tendency to deform the simulation box from its rectangular shape). The relaxation time constants of the thermostat and the barostat were set to 1.0 ps. Lennard-Jones interactions were truncated at 10 Å, and long-range electrostatics were treated with the Particle Mesh Ewald (PME) method65,66 and a k-grid spacing of 0.012 Å−1. We checked that the box vectors and the enthalpies were always well-converged during the equilibration simulations used for the optimizations.

Figure 3. Final conﬁgurations after geometry optimizations using DFT for water molecules on TiO2 surfaces. Panels a and b are side views of the (100) and (110) surfaces, and c and d are the corresponding bottom views. The water molecule is decomposed by exposed Ti atoms but bind noncovalently to surface oxygens.

[(100) surface] and 3b [(110) surface]. Figure 3 shows that exposed Ti surface atoms can decompose the water molecule, leaving the surface hydroxylated (with OH-groups binding to the exposed titaniums). In contrast, the water molecule is kept intact when only oxygens are present at the surface and binds noncovalently to the (100) surface. The distance between the water oxygen and the (100) surface is thus 2.86 Å, but only 1.87 Å for the (110) surface, which is a distance associated with covalent binding. Using the TIP3P model value σOW = 3.16 Å and applying the Berthelot−Lorentz combination rules “backwards” (see “TiO2 Bulk Models” section, but instead given σij and σj, ﬁnd σi), we ﬁnd the approximative values σOT2 = 2.72 Å and σTi5 = 4.303 Å to be used as initial values for the optimization of these parameters. The interpretation of the DFT calculations is that the ﬁrst solvation layer on TiO2 surfaces will range from strongly adsorbed (100 surface) to hydroxylated (110 surface). Such strong adsorption must be accounted for explicitly in classical models by predetermination of the degree of hydroxylation (which is pH-dependent). In our bonded model, we have

■

RESULTS AND DISCUSSION We now proceed to the results of the optimization, and begin with reporting the ﬁndings of the geometry optimizations and 18115

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

Figure 4. (a) The objective function χ2 converges after ∼1000 steps during the optimization of the MA model. The two bottom panels show the target data calculated from the simulations during the optimization; (b) deviation from TiO2 lattice parameters in units of a = 4.592 Å and c = 2.958 Å, and (c) deviation from the TiO2−water adsorption enthalpy.

plateau values were obtained after 1000 iteration steps. The σ’s decreased with ∼5%, which is the most important factor that determines the bulk structure (lattice parameters). Meanwhile, the 14% decrease of ϵTi and the 56% decrease of ϵO is actually a minor contribution to the water adsorption enthalpy compared to the 35% decrease in qTi, because the polar TiO2 material is dominated by the electrostatic interactions. The bulk charge of TiO2 converges to 1.433e, which, interestingly, is in close agreement to the value determined from the experimental electron density data (1.38e) and also to semiempirical calculations based on electronegativity scales. Further, the converged value ϵOT = 0.615 kJ/mol is in agreement with the well-established value ϵO = 0.636 kJ/mol for bulk water, and the converged value ϵTi = 2.181 kJ/mol is close to the value estimated from the London equation (ϵTi = 1.85 kJ/mol). We conclude that all parameters converged to values in agreement what can be expected from physical arguments. Table 3 summarizes the optimization of the MA model, and Figure 4 shows how the target data is reproduced during the optimization. The algorithm continually counterbalances the targets to ﬁnd the best compromise in producing the full set of target data. Table 2 shows that the MA model can be improved to reproduce the adsorption enthalpy to within 15% of the experimental value (from 72%) without additional parameters. The accuracy in the lattice parameters are −1% (a,b) and +4% (c) (from +2% and +6.6%). Optimizing the Bonded TiO2 Model. The most serious approximation in the MA model is the neglect of the covalent contribution to the Ti−O bond, as found in the experimental electron density of rutile.30 The MA model rectiﬁes this error with a large partial charge (qTi = 2.196e) that is counterbalanced by small σ’s and large ϵ’s (to get the correct structure). This cancellation of errors works for modeling bulk phases but fails when simulating interfaces. We now parametrize a TiO2 model that explicitly takes the covalent Ti−O bond into account. We call this the “bonded”

therefore chosen to consider only oxygen-covered TiO 2 surfaces. If surfaces with exposed Ti atoms are simulated, the undercoordinated Ti atoms should be capped with OH-groups prior to the simulation. Improving the Matsui−Akaogi Model. The ﬁrst step was to apply the procedure as-is to the Matsui−Akaogi (MA) model. Our calculations show that the MA model overestimates the TiO2−water enthalpy, Δhads, by 72% compared to experimental data (Table 2). Perhaps this is not too surprising given that the MA model was parametrized to the lattice parameters and elastic constants of the bulk phases of TiO2 (rutile, anatase, and brookite), but not to surface interactions. The MA model does not take the covalent contribution of the Ti−O bond into account, and therefore requires small values of σ and ϵ to obtain the correct TiO2 structure. This imbalance leads to severe errors for surface interactions. However, even though the lattice parameters in the original paper were reported to be accurate to 2% (a,b) and 3% (c), we found the corresponding numbers in our simulations of the MA model to be 2% and 6.6%. This led us to see whether we could improve the MA parameters by systematic optimization of Δhads, a = b, and c in ForceBalance. Figure 4 shows the optimization, which converged in roughly 1000 steps. The noise present when calculating χ2 is inherent to the thermal ﬂuctuations that accompany MD simulations. The algorithm must therefore be willing to accept optimization steps that correspond to small increases in χ2 during the optimization. The noise level is minimized by high-accuracy evaluations of the parameter derivatives that follow from the analytical expressions given by eq 9 and eq 10. The ﬁve force ﬁeld parameters (σ and ϵ for Ti and O, and the TiO2 unit charge q) were varied simultaneously to ﬁnd the optimal set in reproducing the target data. Figure 5 shows the change of the parameters during the optimization. The parameters in the model all decreased, but with varying magnitude (which is also related to the prior); converged 18116

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

Figure 5. Evolution of the force ﬁeld parameters during the optimization of the MA model.

adaptive trust radius with minimum value 0.2. Stage 1 converged within 100 steps, and stage 2 converged within 500 steps. The results of the optimization are summarized in Table 4. Figure 6 shows the objective function during the ﬁrst stage of the optimization of the bonded TiO2 model and Figure 6b,c shows the corresponding target data. The bond length is restricted with a small prior and must stay close to its initial value. The structure is independent of the bond force constant as shown in Table 4. The lattice parameters (and the bulk density) are primarily determined by σTi and σOT, i.e., by the atomic size but excluding 1−2- and 1−3-interactions for bonded neighbors. The two σ’s are similar to the values predicted from the Cambridge Structural Database50 but are adjusted by 5−10% by the optimization procedure. The σ’s converge to their optimal values within 10 iteration steps (Figure 7), after which the objective function is prevented from decreasing further because of the rigid structure imposed by the covalent bonds (visible as oscillations in Figure 6). The lattice parameters a and b versus c can not be simultaneously be made to reproduce their target valus closer than within +1.2% and −3.0%, but the opposite signs cancel the error for the density, which is therefore in excellent agreement with the reference value (−0.65%). Radial distribution functions (RDFs) were calculated to compare the TiO2 structures produced by the diﬀerent models. The RDF between atom types a and b is deﬁned by

Table 3. Parameter Optimization of the Nonbonded TiO2 Modela parameter

unit

α

p0 (MA)

popt (optimized)

%

qTi σTi σOT ϵTi ϵO

e Å Å kJ/mol kJ/mol

1.0 0.5 0.5 0.5 0.5

2.196 1.958 2.875 2.5422 1.3897

1.432 1.8235 2.7486 2.181 0.615

−35 −7 −4 −14 −56

α is the prior; the interval within which the initial parameter is expected to vary. p0 and popt are the initial and optimized values of the parameter. % denotes the deviation between the initial and optimized values.

a

model. We account for surface chemistry by selecting atom types depending on coordination numbers, so that there are 2 atom types for Ti (5- and 6-coordination) and 2 atom types for O (2- and 3-coordination). This approach allows simulations of realistic surface chemistry for nanostructured TiO2 in contact with biomaterials, for the price of a model with more parameters. (The bonded model can not be used to simulate bulk phase transformations). The target data used for the bulk structure (lattice parameters) can be optimized independently from the target data used for the surface interactions (adsorption enthalpy). The optimization of the bonded model then follows the two-stage procedure: 1 Optimizing the lattice parameters (a = b and c) of TiO2 by varying the bulk parameters (Lennard-Jones parameters for the Ti6 and OT3 atom types, and the parameters for the Ti−O bond force constant, kb, and equilibrium distance, b0). We ﬁxed the charge of the TiO2 unit to q = 1.433e in agreement with the optimized MA model (and the experimental value). 2 Optimizing the TiO2−water adsorption enthalpy (Δhads) by varying the surface interaction parameters (LennardJones parameters for the Ti5 and OT2 atom types and the charge qTi5). The bulk parameters are ﬁxed at their optimized values from stage 1. The optimization was performed with the Levenburg− Marquart algorithm as implemented in ForceBalance, with an

gab(r ) =

V NaNb

Na , Nb

∑

δ(r − rijab) (14)

i≠j

rab ij

where the brackets denote the ensemble average, and is the distance between atoms i (type a) and j (type b). The δfunction is 1 if rij lies within the interval r + Δr and 0 otherwise. V is the volume of the simulation box, while Na and Nb are the total number of a and b atoms, respectively. Figure 8 shows the RDFs for the Ti−Ti, Ti−O and O−O pairs in bulk simulations of the MA model, the optimized MA model (nonbonded), and the new bonded model, as compared to the RDF from the 18117

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C Table 4. Parameter Optimization of the Bonded TiO2 Modela parameter

unit

α

p0

popt

%

σTi6 σTi5 σOT3 σOT2 ϵTi6 ϵTi5 ϵO3 ϵO2 qTi6 qTi5 bTi−O kTi−O

Å Å Å Å kJ/mol kJ/mol kJ/mol kJ/mol e e Å kJ/mol/Å2

0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.0 0.5 0.1 1000

4.92 4.92 3.0 2.72 1.815 1.8168 0.636 0.636 1.433 1.8 1.96 2131

5.231 4.6978 3.329 2.6951 1.8168 1.8151 0.6480 0.6226 1.433 1.8365 1.9598 2113

+6.3 −4.5 +11.0 −0.9 +0.1 −0.1 +1.9 −2.1 +0.0 +2.0 −0.0 −0.9

force ﬁeld parameters of the atom types with lower coordination numbers (5-coordinated Ti’s and 2-coordinated O’s). Figure 9 shows the objective function during the second optimization stage; χ2 decreases rapidly the ﬁrst 10 steps, followed by a slow descent until convergence is reached after ∼200 steps. The target data, the adsorption enthalpy Δhads, follows a similar trend as χ2 during the optimization (Figure 9). When the optimization has reached convergence, the calculated target data undulates above the reference value until it is reached at iteration 443 with some help from a ﬂuctuation. This implies that there is a volume in parameter space that corresponds to a range of diﬀerent parameter values, which equally well reproduce the reference target value of the adsorption enthalpy. Figure 10 shows the change of the force ﬁeld parameters during the second stage of the optimization. Keeping in mind that χ2 converges after 250 steps, the parameter with most impact is σTi5 (decreases by −4.5%). The changes in the other parameters are statistically insigniﬁcant, with the partial exception of the 2% increase of q5. Structure and Dynamics of Water at the TiO2 Surface. Figure 11 shows the structure of single adsorbed water molecules on the TiO2 (100) surface obtained with the diﬀerent models. The water molecule penetrates the ﬁrst surface layer because of the small σ-values employed in the nonbonded models (to ensure correct TiO2 bulk structure). The adsorbed water structure obtained with the bonded model is very similar to the DFT structure. The nearest oxygen−oxygen distance with the bonded model is slightly smaller than from the DFT calculation (2.64 vs 2.86 Å), but the structure with hydrogen bonding to the surface oxygen is the same. Water behavior at the TiO2 interface was further investigated for the diﬀerent models by running 500 ns-simulations of fully hydrated TiO2 slabs (∼500 water molecules). The structure of the surface water was analyzed by calculating water density

α is the prior; the interval within which the initial parameter is expected to vary. p0, and popt are the initial and optimized values of the parameter. % denotes the deviation between the initial and optimized values.

a

experimentally determined crystal structure of rutile.The RDF peaks of the simulated models are broadened by thermal ﬂuctuations, which naturally are smallest in the bonded model. The Ti−O distance, rTi−O, is the most important feature and is represented by the innermost peak in the middle panel of Figure 8. rTi−O is overestimated by the MA model, but in good agreement for its optimized counterpart. The best agreement is of course obtained from the bonded model, which restrains the Ti−O distance with a harmonic bond. The overall agreement to the crystal structure RDF is better with the optimized version of the MA model compared to the original version, but the bonded model outperforms the other two. The surface properties of the bonded model are optimized in the second stage of the optimization, and are determined by the

Figure 6. (a) The objective function reaches a minimum in 10 steps during the ﬁrst stage (bulk) of the optimization of the bonded TiO2 model. The lower panels show the target data (in units of each lattice parameter) calculated from the simulations during the ﬁrst stage of the optimization of the bonded model; (b) lattice vector a and (c), lattice vector c. 18118

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

Figure 7. Evolution of the force ﬁeld parameters during the ﬁrst stage of the optimization of the bonded TiO2 model. The varied parameters were the Lennard-Jones parameters for the atom types representing bulk coordination (σOT3, ϵOT3, σTi6 and ϵTi6), the equilibrium Ti−O bond distance (b0) and the Ti−O bond force constant (kb).

Figure 8. Radial distribution functions (RDFs) for pairs of (top) Ti−Ti, (middle) Ti−O, and (bottom) O−O atoms for the bonded (optimized), nonbonded (optimized Matsui−Akaogi) and MA (original Matsui−Akaogi15) models. The dotted lines are the reference RDFs obtained from the crystal structure of rutile.

respect to bulk water. The electric potential was calculated by integrating the charge density twice with respect to the z-direction. Figure 12 shows the density proﬁles and electric potentials for the three TiO2 models. Water was found to strongly adsorb to the TiO2 surface for all three models (MA, nonbonded, and bonded). When the MA force ﬁeld or its optimized nonbonded

proﬁles and electric potentials along the surface normal (the z-direction). The density proﬁles were obtained by dividing the simulation box into N = 1000 slices of equal thickness ΔzLz/N along the z-direction and performing a mass-weighted histogram for the bins. The surface position is taken to coincide with the top layer in the TiO2 slab. The water density proﬁles were shifted to the position of this peak, and then normalized with 18119

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

Figure 9. (a) The objective function χ2 converges after ∼200 steps during the second stage (surface) of the optimization of the bonded TiO2 model. (b) The target data (TiO2−water adsorption enthalpy) calculated from the simulations during the optimization.

Figure 10. Evolution of the force ﬁeld parameters during the second stage of the optimization of the bonded TiO2 model. The varied parameters were the Lennard-Jones parameters for the atom types with surface coordination (σOT2, ϵOT2, σTi5, and ϵTi5p) and the charge of the 5-coordinated Ti atom type (q5, see eq 6).

waters and bulk waters. The penetrating layer is separated from the second solvation layer by hydrogen. The density proﬁles of the MA and the nonbonded models are similar. They show strongly ordered water at the interface (three solvation layers with solid-like properties can easily be distinguished in Figure 12) with a broad shoulder that extend almost 10 Å into the bulk. The force ﬁeld parameters of the nonbonded model

counterpart is used to model TiO2, water molecules penetrate the ﬁrst surface layer (shown by the ﬁrst peak in the density proﬁles in Figure 12). This is possible because of the small σ-values (atomic sizes) of Ti and O in these two models. Penetrating water molecules stay bonded to the surface for the entire 500 ns simulation and will be referred to as coordinated water molecules. No exchange occurs between coordinated 18120

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

Figure 11. Conﬁgurations of adsorbed water molecules on the TiO2 (100) surface obtained with the (a) MA model, (b) nonbonded model, and (c) bonded model. The conﬁguration obtained from geometry optimizations using DFT is shown with transparent rendering. The water molecule penetrates the ﬁrst surface layer when the nonbonded models (a and b) are used, while the water structure obtained with the bonded model (c) is similar to the DFT structure (with hydrogen bonding to the surface oxygens).

Figure 12. Panels of each row corresponds to simulations with the MA (top), nonbonded (middle) and bonded (bottom) TiO2 models. The panels in the left column shows the mass density of the water molecule and the water oxygens, in the normal direction to the surface. The mass density has been normalized with respect to the water bulk density. z = 0 corresponds to the top layer of the TiO2 slab. The innermost peaks of the MA and nonbonded models are coordinated waters. The panels in the right column right column shows the electrical potential (solid lines) and the water mass density (dash-dotted lines) in the normal direction to the surface. The arrows mark the values used to calculate the ζ-potential. The distance to the surface corresponds to the ﬁrst liquid water layer.

that the electric potential oscillates due to the ordered water structure near the surface. The potential diﬀerence between the solvent and the stationary layer of ﬂuid attached to the TiO2 surface (the “slipping plane”) is known as the ζ-potential and can be measured experimentally, e.g., with electrophoresis.67 The ζ-potential for TiO2 decreases with increasing pH until leveling out around −50 mV at neutral pH 7 (see e.g., Figure 7 in ref 39). Deﬁning the slipping plane (which separates surfacebound ﬂuid from bulk) unambiguously from simulations is nontrivial. Interestingly, we ﬁnd a minimum in the electrical potential near the surface that coincides with the forth (nonbonded models) or third (bonded) solvation layer (see Figure 12). Using this minimum to deﬁne the slipping plane, we ﬁnd the ζ-potentials −290 mV (MA model), −145 mV

are less electrostatic (polar) than those of the original MA model, but the same functional form of the energy expression (no covalent contribution to the Ti−O bond which is compensated by small atomic σ’s). Therefore, the density proﬁles of the original and the optimized MA models are similar. In contrast, the density proﬁle of the bonded TiO2 model is different. There is a single peak of adsorbed waters around ∼2 Å from the surface. This water layer is solid-like (ice) but not irreversibly bound to the surface; there is exchange between surface and bulk waters. A second solvation layer is shown in the density proﬁle but the magnitude of the ordering is far less compared to the MA and nonbonded models. There is no oriental ordering in bulk water (z → ∞) and the electric potential averages to zero. However, Figure 12 shows 18121

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C

Figure 13. C(t) is the contact autocorrelation function for the bonded TiO2 model and measures how long water molecules stay adsorbed to the TiO2 surface. The solid line is an exponential ﬁt to simulation data (circles) with residence time τ = 31 ns.

(nonbonded model) and −70 mV (bonded model). The ﬁrst two values are overestimates in comparison to the experimental value (−50 mV at neutral pH), while the bonded model is in favorable agreement with experimental data. Note that this negative potential originates in speciﬁc orientations of water dipoles at the surface, and not in explicit surface charge due to protonation. Further validation of TiO2/water interactions can be determined by comparison to heats of immersions as determined in experiments: ΔHimm = (Hdry + H water − H wet)/A

however, provides realistic values for all the TiO2−water interactions studied here. The interfacial dynamics in the bonded model was studied by monitoring how long water molecules stay adsorbed to the surface. Adsorption is determined by the “contact function” ni(t), deﬁned by ⎧ 1 when molecule i is adsorbed ni(t ) = ⎨ ⎩ 0 otherwise

(16)

for (water) molecule i. Here, the molecule is adsorbed when the water oxygen is within 3 Å of the TiO2 surface, which is the distance to the ﬁrst minimum in the radial distribution function of TiO2 and OW (water oxygen). We calculated the contact autocorrelation function (CAF) according to

(15)

where Hwet is the enthalpy of a TiO2 slab solvated by N water molecules, Hwater is the enthalpy of N bulk waters, and Hdry is the enthalpy of the TiO2 slab without water. A is twice the lateral area of the surface slab (because of the two interfaces present in simulations with periodic boundary conditions). The experimental range of ΔHimm is 0.2−0.6 J/m2, and is strongly dependent on pretreatment and whether adsorbed waters are present on the dry surfaces.68 The trend is that the heats of immersion increase as the amount of waters present on the dry surfaces decrease, so we should compare with the highest values (since our “dry” surfaces are completely dry). We have calculated ΔHimm to 5.78, 5.06, and 0.79 J/m2 for the MA, nonbonded and bonded models, respectively (using completely dry surfaces). The nonbonded models overestimate the immersion enthalpies by an order of magnitude, while the bonded model is in good agreement with the experimental values. The data (water densities, ζ-potentials, and heats of immersion) shows that the optimization of the nonbonded model improves the target data (the adsorption enthalpy) but fails to reproduce other water interactions. The underlying physical reason is the neglect of the covalent contribution to the Ti−O bond. The ﬂawed ionocovalent balance in nonbonded TiO2 models is counterpoised with small σ-values, which yields a correct TiO2 structure, but too strong interactions when the water interface is present. The bonded model,

C(t ) =

1 N

N

∑ Ci(t ) = i=1

1 N

N

∑ i=1

⟨ni(t )ni(0)⟩ ⟨ni(0)ni(0)⟩

(17)

where N is the total number of water molecules. C(t) is the probability that water molecule i is adsorbed to the surface at time t = t given that it was adsorbed at time t = 0. C(t) starts out at 1 (full correlation) and decays to zero (no correlation) as t → ∞. Note also that C(t) ≥ 0, ∀t by deﬁnition. The CAF can be used to deﬁne a “residence time”, tr =

∫0

∞

dt C(t )

(18)

which is a typical time period a water molecule stays adsorbed to the surface. The residence time coincides with the relaxation time, tr = τ, when the CAF decays exponentially, i.e., C(t ) = e−t/ τ . tr then corresponds to the time when C(t) has decayed to e−1 ≈ 0.3679 of its original value. Figure 13 shows the CAF for the bonded model, which is excellently ﬁtted to exponential decay with relaxation time τ = 31 ns. Exponential decay is characteristic of a two-state model (the molecule can be bound to the surface or free in the 18122

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

The Journal of Physical Chemistry C

■

bulk) with an equilibrium rate between the states that is determined by the relaxation time τ. In a two-state model, the 31 ns relaxation time corresponds to a free energy diﬀerence (adsorption free energy) of water of 30 kJ/mol.69,70

CONCLUSIONS This work has developed accurate molecular models of TiO2 for simulations of the nano−bio interface, using a method for force ﬁeld optimization that accepts arbitrary experimental target data such as structural (lattice parameters) and thermodynamical (adsorption enthalpy) material properties. Properly used, the method furnishes the path to reproducible force ﬁelds developed with automated procedures. The objective function to be minimized can be any combination of experimental or ab initio theoretical data. The method was shown to improve the widely used Matsui− Akaogi force ﬁeld for TiO2 using no additional ﬁt parameters. In addition, a TiO2 model that includes the covalent contribution to the Ti−O bond was developed for surface simulations. The model was based on geometries obtained from periodic DFT calculations and optimized to experimental data. Simulations of fully hydrated TiO2 slabs were used to compare structure and dynamics of interfacial waters for three TiO2 models. The nonbonded models (MA and its optimized counterpart) were shown to strongly order the solvent and lead to “coordinated” waters penetrating into the surface. The bonded TiO2 model has weaker interactions leading to ζ-potential and heat of immersion in agreement with experimental data, and an adsorbed solvation layer that is exchanged with the bulk on the nanosecond-time scale. The TiO2−water structures obtained with the bonded model are in agreement with the optimized geometries from the DFT calculations (Figure 3), with no “coordinated” waters on the TiO2 (100) surface. The bonded TiO2 model is therefore suitable for simulations of interfaces between nanosized TiO2 and biomolecules in explicit solvent. The nonbonded model can be used in phase transition simulations of TiO2 systems, while potential applications for the bonded model include protein adsorption on TiO2 surfaces and corona formation on TiO2 nanoparticles. ASSOCIATED CONTENT

S Supporting Information *

The ﬁnal optimized force ﬁeld parameters for the bonded model and details on performing optimizations in ForceBalance based on thermodynamical data. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b02669.

■

REFERENCES

(1) Ruoslahti, E.; Bhatia, S. N.; Sailor, M. J. Targeting of Drugs and Nanoparticles to Tumors. J. Cell Biol. 2010, 188, 759−768. (2) Kasemo, B. Biological Surface Science. Surf. Sci. 2002, 500, 656− 677. (3) Mahmoudi, M.; Lynch, I.; Ejtehadi, M. R.; Monopoli, M. P.; Bombelli, F. B.; Laurent, S. Protein-Nanoparticle Interactions: Opportunities and Challenges. Chem. Rev. 2011, 111, 5610−5637. (4) Monopoli, M. P.; Åberg, C.; Salvati, A.; Dawson, K. A. Biomolecular Coronas Provide the Biological Identity of Nanosized Materials. Nat. Nanotechnol. 2012, 7, 779−786. (5) Nel, A. E.; Madler, L.; Velegol, D.; Xia, T.; Hoek, E. M. V.; Somasundaran, P.; Klaessig, F.; Castranova, V.; Thompson, M. Understanding Biophysicochemical Interactions at the Nano-Bio Interface. Nat. Mater. 2009, 8, 543−557. (6) Bhattacharya, J.; Choudhuri, U.; Siwach, O.; Sen, P.; Dasgupta, A. K. Interaction of Hemoglobin and Copper Nanoparticles: Implications in Hemoglobinopathy. Nanomedicine 2006, 2, 191−199. (7) Jans, H.; Liu, X.; Austin, L.; Maes, G.; Huo, Q. Dynamic Light Scattering as a Powerful Tool for Gold Nanoparticle Bioconjugation and Biomolecular Binding Studies. Anal. Chem. 2009, 81, 9425−9432. (8) Monopoli, M. P.; Walczyk, D.; Campbell, A.; Elia, G.; Lynch, I.; Baldelli Bombelli, F.; Dawson, K. A. Physical-Chemical Aspects of Protein Corona: Relevance to in Vitro and in Vivo Biological Impacts of Nanoparticles. J. Am. Chem. Soc. 2011, 133, 2525−2534. (9) Lévy, R.; Thanh, N. T. K.; Doty, R. C.; Hussain, I.; Nichols, R. J.; Schiffrin, D. J.; Brust, M.; Fernig, D. G. Rational and Combinatorial Design of Peptide Capping Ligands for Gold Nanoparticles. J. Am. Chem. Soc. 2004, 126, 10076−10084. (10) Cedervall, T.; Lynch, I.; Lindman, S.; BerggÅrd, T.; Thulin, E.; Nilsson, H.; Dawson, K. A.; Linse, S. Understanding the NanoparticleProtein Corona Using Methods to Quantify Exchange Rates and Affinities of Proteins for Nanoparticles. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 2050−2055. (11) Lundqvist, M.; Sethson, I.; Jonsson, B.-H. Protein Adsorption onto Silica Nanoparticles: Conformational Changes Depend on the Particles’ Curvature and the Protein Stability. Langmuir 2004, 20, 10639−10647. (12) Shang, L.; Wang, Y.; Jiang, J.; Dong, S. pH-Dependent Protein Conformational Changes in Albumin:Gold Nanoparticle Bioconjugates: A Spectroscopic Study. Langmuir 2007, 23, 2714−2721. (13) Brandt, E. G., Lyubartsev, A. J. Phys. Chem. C, submitted 2015.10.1021/acs.jpcc.5b02670 (14) Weir, A.; Westerhoff, P.; Fabricius, L.; Hristovski, K.; von Goetz, N. Titanium Dioxide Nanoparticles in Food and Personal Care Products. Environ. Sci. Technol. 2012, 46, 2242−2250. (15) Matsui, M.; Akaogi, M. Molecular Dynamics Simulation of the Structural and Physical Properties of the Four Polymorphs of TiO2. Mol. Simul. 1991, 6, 239−244. (16) Bandura, A. V.; Kubicki, J. D. Derivation of Force Field Parameters for TiO2-H2O Systems from ab Initio Calculations. J. Phys. Chem. B 2003, 107, 11072−11081. (17) Prědota, M.; Bandura, A. V.; Cummings, P. T.; Kubicki, J. D.; Wesolowski, D. J.; Chialvo, A. A.; Machesky, M. L. Electric Double Layer at the Rutile (110) Surface. 1. Structure of Surfaces and Interfacial Water from Molecular Dynamics by Use of ab Initio Potentials. J. Phys. Chem. B 2004, 108, 12049−12060. (18) Alimohammadi, M.; Fichthorn, K. A. A Force Field for the Interaction of Water with TiO2 Surfaces. J. Phys. Chem. C 2011, 115, 24206−24214. (19) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (20) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616.

■

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] Notes

The authors declare no competing ﬁnancial interest.

■

ACKNOWLEDGMENTS This research was supported by the European Commission MembraneNanoPart FP7 project and by the Swedish Research Council (Vetenskapsradet). The authors are thankful to the Swedish National Infrastructure for Computing (SNIC) for granting computer facilities. 18123

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C (21) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (22) Köppen, S.; Langel, W. Simulation of the Interface of (100) Rutile with Aqueous Ionic Solution. Surf. Sci. 2006, 600, 2040−2050. (23) Kim, S.-Y.; Kumar, N.; Persson, P.; Sofo, J.; van Duin, A. C. T.; Kubicki, J. D. Development of a ReaxFF Reactive Force Field for Titanium Dioxide/Water Systems. Langmuir 2013, 29, 7838−7846. (24) Schneider, J.; Ciacchi, L. C. A Classical Potential to Model the Adsorption of Biological Molecules on Oxidized Titanium Surfaces. J. Chem. Theory Comput. 2011, 7, 473−484. (25) Wang, L.-P.; Martinez, T. J.; Pande, V. S. Building Force Fields: An Automatic, Systematic, and Reproducible Approach. J. Phys. Chem. Lett. 2014, 5, 1885−1891. (26) Wang, L.-P.; Head-Gordon, T.; Ponder, J. W.; Ren, P.; Chodera, J. D.; Eastman, P. K.; Martinez, T. J.; Pande, V. S. Systematic Improvement of a Classical Molecular Model of Water. J. Phys. Chem. B 2013, 117, 9956−9972. (27) Materials Studio, version 7.0; Accelrys Software, Inc.: San Diego, CA, 2001−2011. (28) Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J.; Probert, M. J.; Refson, K.; Payne, M. C. First Principles Methods Using CASTEP. Z. Kristallogr. 2005, 220, 567−570. (29) Perdew, J. P., Ziesche, P., Eschrig, H. Electronic Structure of Solids ’91; Akademie Verlag: Berlin, 1991; Vol. 11. (30) Jiang, B.; Zuo, J. M.; Jiang, N.; O’Keeffe, M.; Spence, J. C. H. Charge Density and Chemical Bonding in Rutile, TiO2. Acta Crystallogr., Sect. A: Found. Crystallogr. 2003, 59, 341−350. (31) Duffy, J. A. Ionic-Covalent Character of Metal and Nonmetal Oxides. J. Phys. Chem. A 2006, 110, 13245−13248. (32) Lebouteiller, A.; Courtine, P. Improvement of a Bulk Optical Basicity Table for Oxidic Systems. J. Solid State Chem. 1998, 137, 94− 103. (33) Traylor, J. G.; Smith, H. G.; Nicklow, R. M.; Wilkinson, M. K. Lattice Dynamics of Rutile. Phys. Rev. B 1971, 3, 3457−3472. (34) Heinz, H.; Lin, T.-J.; Kishore Mishra, R.; Emami, F. S. Thermodynamically Consistent Force Fields for the Assembly of Inorganic, Organic, and Biological Nanostructures: The INTERFACE Force Field. Langmuir 2013, 29, 1754−1765. (35) Healy, K. E.; Ducheyne, P. Hydration and Preferential Molecular Adsorption on Titanium in Vitro. Biomaterials 1992, 13, 553−561. (36) Lausmaa, J.; Kasemo, B.; Mattsson, H. Surface Spectroscopic Characterization of Titanium Implant Materials. Appl. Surf. Sci. 1990, 44, 133−146. (37) Roddick-Lanzilotta, A. D.; Connor, P. A.; McQuillan, A. J. An In Situ Infrared Spectroscopic Study of the Adsorption of Lysine to TiO2 from an Aqueous Solution. Langmuir 1998, 14, 6479−6484. (38) Jones, F. H. Teeth and Bones: Applications of Surface Science to Dental Materials and Related Biomaterials. Surf. Sci. Rep. 2001, 42, 75−205. (39) Suttiponparnit, K.; Jiang, J.; Sahu, M.; Suvachittanont, S.; Charinpanitkul, T.; Biswas, P. Role of Surface Area, Primary Particle Size, and Crystal Phase on Titanium Dioxide Nanoparticle Dispersion Properties. Nanoscale Res. Lett. 2011, 6, 27−33. (40) Fitts, J. P.; Machesky, M. L.; Wesolowski, D. J.; Shang, X.; Kubicki, J. D.; Flynn, G. W.; Heinz, T. F.; Eisenthal, K. B. SecondHarmonic Generation and Theoretical Studies of Protonation at the Water/α-TiO2 (1 1 0) Interface. Chem. Phys. Lett. 2005, 411, 399− 403. (41) Machesky, M. L.; Předota, M.; Wesolowski, D. J.; Vlcek, L.; Cummings, P. T.; Rosenqvist, J.; Ridley, M. K.; Kubicki, J. D.; Bandura, A. V.; Kumar, N.; et al. Surface Protonation at the Rutile (110) Interface: Explicit Incorporation of Solvation Structure within the Refined MUSIC Model Framework. Langmuir 2008, 24, 12331− 12339. (42) Israelachvili, J. N. Intermolecular and Surface Forces, 3rd ed.; Academic Press: Waltham, MA, 2011.

(43) Cheng, J.; Sprik, M. Acidity of the Aqueous Rutile TiO2 (110) Surface from Density Functional Theory Based Molecular Dynamics. J. Chem. Theory Comput. 2010, 6, 880−889. (44) Casarin, M.; Maccato, C.; Vittadini, A. Molecular Chemisorption on TiO2 (110): A Local Point of View. J. Phys. Chem. B 1998, 102, 10745−10752. (45) Heinz, H.; Suter, U. W. Atomic Charges for Classical Simulations of Polar Systems. J. Phys. Chem. B 2004, 108, 18341− 18352. (46) Hirshfeld, F. L. Bonded-Atom Fragments for Describing Molecular Charge Densities. Theor. Chim. Acta 1977, 44, 129−138. (47) Coppens, P. X-Ray Charge Densities and Chemical Bonding; International Union of Crystallography Texts on Crystallography; Oxford Science Publications: New York, 1997; Vol. 4. (48) Bondi, A. van der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441−451. (49) Batsanov, S. Van der Waals Radii of Elements. Inorg. Mater. 2001, 37, 871−885. (50) Alvarez, S. A Cartography of the van der Waals Territories. Dalton Trans. 2013, 42, 8617−8636. (51) Heinz, H.; Vaia, R. A.; Farmer, B. L.; Naik, R. R. Accurate Simulation of Surfaces and Interfaces of Face-Centered Cubic Metals Using 12−6 and 9−6 Lennard-Jones Potentials. J. Phys. Chem. C 2008, 112, 17281−17290. (52) Lide, D. R., Ed. CRC Handbook of Chemistry and Physics, 86th ed.; CRC Press: Boca Raton, FL, 2005. (53) Grant, F. A. Properties of Rutile (Titanium Dioxide). Rev. Mod. Phys. 1959, 31, 646−674. (54) Diebold, U. The Surface Science of Titanium Dioxide. Surf. Sci. Rep. 2003, 48, 53−229. (55) Downs, R. T. Program and Abstracts, 19th General Meeting of the International Mineralogical Association, Kobe, Japan, 2006. (56) Lyubartsev, A.; Mirzoev, A.; Chen, L.-J.; Laaksonen, A. Systematic Coarse-Graining of Molecular Models by the Newton Inversion Method. Faraday Discuss. 2010, 144, 43−56. (57) Swope, R. J.; Smyth, J. R.; Larson, A. C. H in Rutile-Type Compounds; I, Single-Crystal Neutron and X-Ray Diffraction Study of H in Rutile. Am. Mineral. 1995, 80, 448−453. (58) Bolis, V.; Busco, C.; Ciarletta, M.; Distasi, C.; Erriquez, J.; Fenoglio, I.; Livraghi, S.; Morel, S. Hydrophilic/Hydrophobic Features of TiO2 Nanoparticles as a Function of Crystal Phase, Surface Area and Coating, in Relation to Their Potential Toxicity in Peripheral Nervous System. J. Colloid Interface Sci. 2012, 369, 28−39. (59) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (60) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (61) Momma, K.; Izumi, F. VESTA3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. J. Appl. Crystallogr. 2011, 44, 1272−1276. (62) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (63) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to An External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (64) Ming, L.-C.; Manghnani, M. H. Isothermal Compression of TiO2 (Rutile) Under Hydrostatic Pressure to 106 Kbar. J. Geophys. Res. 1979, 84, 4777−4779. (65) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (66) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N × log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. 18124

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125

Article

The Journal of Physical Chemistry C (67) Delgado, A. V.; González-Caballero, F.; Hunter, R. J.; Koopal, L. K.; Lyklema, J. Measurement and Interpretation of Electrokinetic Phenomena. J. Colloid Interface Sci. 2007, 309, 194−224. (68) Morimoto, T.; Nagao, M.; Omori, T. Heat of Immersion of Titanium Dioxide in Water. I. The Effect of the Hydration Treatment of Titanium Dioxide. Bull. Chem. Soc. Jpn. 1969, 42, 943−946. (69) Markovitch, O.; Agmon, N. Reversible Geminate Recombination of Hydrogen-Bonded Water Molecule Pair. J. Chem. Phys. 2008, 129, 084505. (70) van der Spoel, D.; van Maaren, P. J.; Larsson, P.; Tîmneanu, N. Thermodynamics of Hydrogen Bonding in Hydrophilic and Hydrophobic Media. J. Phys. Chem. B 2006, 110, 4393−4398.

18125

DOI: 10.1021/acs.jpcc.5b02669 J. Phys. Chem. C 2015, 119, 18110−18125