Benchmarks and Dielectric Constants for Reparametrized OPLS and

Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

Article

Benchmarks and Dielectric Constants for Reparametrized OPLS and Polarizable Force Field Models of Chlorinated Hydrocarbons Zhu Liu, Jakob Timmermann, Karsten Reuter, and Christoph Scheurer J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b06709 • Publication Date (Web): 23 Aug 2017 Downloaded from http://pubs.acs.org on August 23, 2017

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

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

Page 1 of 31

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

The Journal of Physical Chemistry

Benchmarks and Dielectric Constants for Reparametrized OPLS and Polarizable Force Field Models of Chlorinated Hydrocarbons Zhu Liu, Jakob Timmermann, Karsten Reuter, and Christoph Scheurer*

Chair for Theoretical Chemistry and Catalysis Research Center, Technische Universität München, Lichtenbergstr. 4, 85747 Garching, Germany

ABSTRACT The knowledge of dielectric response properties of the environment is of paramount importance in many theoretical embedding methods and studies of solutes and of catalytic sites and processes in condensed phases. In particular, the realistic embedding of active sites into solid/liquid and liquid/liquid interfaces is a crucial point in the context of modelling energy conversion (e.g. electrochemical, photochemical, power-to-X) processes. Recently, the finding that the dielectric permeability of liquids near solid/liquid interfaces is far from being constant but deviates strongly from the bulk value within several nanometers from the interface has raised the interest in a more fundamental understanding of the response properties near interfaces. As these questions are hard to study experimentally, reliable theoretical models are required. Here we describe a careful first-principles based reparametrization of nonpolarizable molecular mechanics force fields for a class of technological relevant organic chlorinated hydrocarbon solvents which are immiscible with water. For the solvent 1,2-dichloroethane (1,2-DCE) we also present a new polarizable force field based on the Drude oscillator model. Its parametrization needs particular attention to avoid unphysical couplings between the internal torsional degree of freedom and the Drude oscillators, which could severely skew the response properties. The performance of this new set of force fields is critically assessed based on a comprehensive molecular dynamics study.

ACS Paragon Plus Environment

1


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

Page 2 of 31

1. INTRODUCTION Chlorinated hydrocarbons (CHCs) are a very large and diverse group of organic molecules characterized by their hydrocarbon skeleton with at least one chlorine atom chemically bonded to it. Notwithstanding that some of its members environmental impact has to be rated critically, this class of compounds has seen widespread use e.g. as solvents or chemical precursors in numerous application fields that are of great economic and technological relevance.1-2 Many of the CHCs with one or two carbon and a varying number of chlorine atoms are rather nonpolar solvents with low boiling points and dielectric constants that are immiscible with water. The formation of a liquid/liquid interface between water and a CHC is particularly interesting for those cases where the CHC phase can be “polarized” substantially in the vicinity of the interface by the presence of the water phase. These cases have hence drawn attention also in theoretical studies.3-5 In particular, 1,2-dichloroethane (1,2-DCE) constitutes an interesting prototype model molecule, owing to the considerable difference in the dipole moments of its gauche and trans conformers.5-11 Both conformers occur in a dynamic equilibrium in the liquid phase of 1,2-DCE and the trans conformer’s dipole moment vanishes due to its inversion symmetry. The detailed microscopic understanding of phase transfer or catalytic and (photo-)electrochemical processes at complexes embedded in such liquid/liquid (but also solid/liquid) interfaces has to rely heavily on simulations12-15 for which the dielectric response of the involved phases is one of the key properties.16-20 The literature on water models and their dielectric response behavior is understandably vast, in particular compared to the CHCs. Reliable modeling of liquid/liquid interfaces requires an equally sound understanding of both involved phases, though. For mechanical and thermal properties of the pure CHC liquid phases, molecular dynamics (MD) simulations have proven to be a reliable tool to understand the behavior of these organic liquids at the molecular level and to link the microscopic behavior statistically to macroscopic properties.4-6 The accurate calculation of thermodynamic and dynamic properties of organic molecular fluids via MD simulations depends crucially on the reliability of the effective molecular mechanics force field to describe all relevant


2

Page 3 of 31

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


intramolecular and intermolecular interactions. This will be even more important for the simulation of liquid/liquid interfaces of CHCs in contact with a liquid water phase. Over the past decades, several force fields such as OPLS-AA (Optimized Potentials for Liquid Simulations – All Atom)21, GAFF (Generalized Amber force field)22-23, and CHARMM24 have been developed for solutes and solvents to model solution systems by MD. E.g. the group of Jorgensen8, 21, 25-27 showed that the OPLS-AA force field is generally well suited for the simulation of a wide variety of nonaqueous solvents. Recently, Caleman et al.23, 28 presented an extensive MD-based study to provide a force field benchmark for organic liquids comparing the OPLS-AA and GAFF models. Target properties like density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and the bulk dielectric constant of a set of organic liquids were compared with experimental data. Most properties were in reasonably good agreement, with the notable exception of the dielectric constant which exhibited rather large deviations from the experimental values. The systematic procedure applied by Frank to parametrize force fields for molecular fluids can achieve improved dielectric constants, but only fully rigid MD models were included.7 While rigid models are computationally more efficient, only flexible models will be able to capture the intramolecular dynamics faithfully, which is necessary to simulate the physical and dielectric properties e.g. of the 1,2-DCE molecule described above. In the current work we will thus focus exclusively on flexible all-atom force fields. Due to its simplicity and computational efficiency on the one hand and its success in reproducing numerous thermodynamic, structural, and dynamical properties of the CHCs on the other hand, we choose the OPLS-AA force field as a starting point. Given our current knowledge about the CHCs, however, the original OPLS-AA force field is known to have difficulties predicting sufficiently accurate values of bulk dielectric constants.23 The systematic improvement of the dielectric constant computed from MD simulations with the reparametrized CHC models is thus one of the key points of this paper. Nonpolarizable (NP) models, like OPLS-AA, are still commonly used due to the relatively simple potential energy functions, low computational costs, and their


3


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

Page 4 of 31

remarkable success in modeling many molecular systems. In general, for improving on an efficient force field, we should aim for strategies to increase its accuracy without increasing the number of interaction sites.7, 29-30 For most of the CHCs we thus limit ourselves to the NP OPLS-AA type of model, where this seems sufficiently accurate. For water, the importance of accounting for molecular polarizability in modeling aqueous phenomena has been demonstrated long ago.31 Incorporating electronic polarization is thus seen as a major requirement for next generation force fields.3 For 1,2-DCE we also present a polarizable force field based on the classical Drude oscillator (DO) model32-35 in order to further improve the prediction of the physical properties in comparison to the NP model. The DO model is an efficient tool to account for electronic polarizability without incurring significantly increased computational costs.32-33 We adopt the approach by Lamoureux and Roux32 and introduce polarizable carbon atom sites. The self-consistent field (SCF) treatment of the polarizability is avoided by treating the DO as an additional classical degree of freedom in the extended Lagrangian mechanics formalism.36-37 For the DO model as well as for the NP force field reparametrization, we focus in particular on improved DCE torsional potential energy profiles as compared to first-principles calculations and partial charges based on the RESP38-39 fit approach for an improved description of the molecular dipole moments. A careful and balanced treatment of these properties forms the basis for improving upon earlier CHC models in the description of the dielectric response. This paper is organized as follows. In Section 2, the method to derive the reparametrization is briefly described, followed by the details of the MD simulations. In Section 3, the analysis procedures used to compute the target properties of the improved models are outlined. Section 4 presents the resulting properties for the improved CHC models and including comparisons with the available experimental data40-42 as well as the work of the Caleman et al.23, 28 on the OPLS-AA force field. A conclusion is given in Section 5.

2. METHODS The OPLS-AA force field studied by Caleman et al.23, 28 was used as a starting point


4

Page 5 of 31

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


for our NP reparametrizations. The standard bond stretching, angle bending and nonbonded Lennard-Jones parameters were initially assigned from the OPLS-AA parameter set. The present work then focused on the improvement of the torsional angle coefficients (where applicable) and atomic partial charges within the flexible NP force field. The objective was to produce simple and efficient OPLS-AA type CHC models for chloroform (CHCl3), dichloromethane (CH2Cl2), 1,1-dichloroethane (1,1-DCE), and 1,2-dichloroethane (1,2-DCE-NP) that yield improved bulk dielectric constants while at least retaining or possibly improving the computed values for all other relevant physical properties of the solvents. The same approach was then followed to parametrize a polarizable model for 1,2-DCE (1,2-DCE-DO, “DO” for Drude oscillator), now including polarizable carbon atoms to further improve molecular bulk properties. As target properties of the liquids besides the bulk dielectric constant, we focused on equilibrium values for the density, heat capacities, isothermal compressibility, and the volumetric expansion coefficient. Properties of the CHC liquid-vapor interfaces such as the heat of vaporization and surface tension were evaluated in addition.

2.1 Nonpolarizable reparametrization 2.1.1 Atomic charges The atomic partial charges for the different interaction sites were fitted using the RESP38-39 method. The molecular electrostatic reference potential was computed using GAUSSIAN 03 with the 6-311G* basis set43, in order to be in line with the original parametrization of the OPLS force field21. For 1,2-DCE-NP, the atomic charge distributions were slightly refitted to achieve better agreement with the respective dipole moment profile as a function of the Cl-C-C-Cl torsion angle. For this purpose, the negative charge of the 1,2-DCE-NP carbon was slightly increased (-0.0075) while the charge of the chlorine was slightly decreased (+0.0075), which result in the cis conformer model of the 1,2-DCE-NP having a lower over-all dipole moment. The same charge modification was applied to the two different carbon atoms of the 1,1-DCE model but without further fitting, increasing its over-all dipole


5


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

Page 6 of 31

moment. For CHCl3, the carbon charge was fixed at zero to avoid overfitting, resulting in non-physical partial charges. Additionally, the models for CHCl3 and CH2Cl2 with the plain RESP38-39 charges did not represent the target properties well enough, requiring slight readjustments of the atomic partial charges based on short MD simulations to achieve better agreement with the experimental properties. The final atomic charges for all interaction sites of the CHCs are collected in Table 1.

2.1.2 Torsion potential parameters For 1,1-DCE and 1,2-DCE, the torsional potential energy profiles depending on dihedral angles, φi, are described by cosine series of the form27 1 1 1 1 Etorsion = ∑ [ V1,i (1 + cos ϕi ) + V2,i (1 − cos 2ϕi ) + V3,i (1 + cos 3ϕi ) + V4,i (1 − cos 4ϕi )] . 2 2 2 2 i The torsion parameters Vk,i were fitted to the first-principles energy profiles using the differential evolution algorithm from the Python package INSPYRED 1.044. The torsional cut through the first-principles reference potential energy surface was evaluated at 19 different 1,2-DCE conformer geometries with GAUSSIAN at the B3LYP-aug-cc-pVDZ45 level, scanning the dihedral angle from 0 to 180. To be able to fit the torsional energy cut for non-polarizable as well as for polarizable models (1,2-DCE-DO) on equal footing, a common protocol was followed. For each of the target conformer geometries a single molecule MD simulation was conducted, progressively cooling down the system while simultaneously restraining the dihedral angle and eventually extracting an effective torsional potential energy value for the respective conformer. The conformer geometries were further weighted properly to keep the focus on the two main conformers (gauche and trans) and the relative height of the lower one of the two rotational energy barriers. Out of an ensemble of parameter sets generated by this protocol with comparable fit quality for the energy profile the one set was chosen that reproduces the experimentally observed trans/gauche distribution best within an MD simulation. The resulting numerical coefficients for this torsion potential energy are listed in the Table 2.


6

Page 7 of 31

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


Table 1. Partial atomic charges q for the interaction sites of the CHC molecules (* value is restrained by construction, see text).

Molecule

Interaction site

q (e)

CHCl3

C

0.0000*

H

0.0700

Cl

-0.2100

C

-0.2140

H

0.2075

Cl

-0.1005

C(H3C-)

-0.1615

H(CH3-)

0.0950

C(Cl2HC-)

-0.0525

H(Cl2CH-)

0.1920

Cl

-0.1315

C

-0.0515

H

0.1160

Cl

-0.1805

DC

2.1911

H

0.1160

Cl

-0.1805

DP

-2.2426

CH2Cl2

1,1-DCE

1,2-DCE-NP

1,2-DCE-DO

Table 2. Coefficients of the torsion potential energy for 1,1-DCE and 1,2-DCE.

Molecule

Dihedral angle

V1

V2

V3

V4

1,1-DCE

H-C-C-H

1.5000

0.2960

0.1910

0.1080

H-C-C-Cl

0.2300

0.9120

0.3590

0.0000

H-C-C-H

1.5000

0.2960

0.1910

0.1080

H-C-C-Cl

0.2300

0.9120

0.3590

0.0000

Cl-C-C-Cl

1.4490

0.4960

0.1270

0.0000

1,2-DCE-NP


7


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

1,2-DCE-DO

Page 8 of 31

H-C-C-H

1.5000

0.0000

0.1603

1.2498

H-C-C-Cl

0.0696

0.6973

0.2454

0.6100

Cl-C-C-Cl

1.5000

0.2546

1.0054

0.0614

2.2 Polarizable reparametrization As mentioned above, NP models are not capable of fully capturing the molecular response to external fields due to their lack of a polarizable electron density. One way to take this contribution into account approximately is given by the Drude oscillator model.34-35 This model accounts for the induced polarization by introducing an almost massless Drude particle (DP) that is harmonically attached to a massive core (DC). The sum of the masses of the DC-DP pair should equal the mass of the modelled atom. Analogously, DC and DP carry different partial charges summing up to the charge of the target atom. The positions of these particles are relaxed into their local minimum energy position in a self-consistent field (SCF) manner, which has to be solved at each MD time step, making the computational cost increase noticably.32 What makes the DO model computationally tractable in MD simulation is the treatment of the Drude oscillators as classical dynamic variables in the context of the extended Lagrangian methods.36-37 This is achieved by allocating a small mass mD to the DP from the initial atom and applying separate thermostats to the relative motion of the DPs around their DCs and the center of mass of the DC-DP pairs.36-37 In the current model with a polarizable carbon, mD is chosen to be small, mD = 0.4 g/mol, and the force constant kD of the harmonic restoring potential is set to be quite large, kD = 4180kJ/mol/Å2, for all Drude oscillators in order to maintain the kinetic decoupling of the DC-DP pair motion from the other degrees of freedom. In addition, two separate Langevin thermostats are applied to thermalize the reduced degrees of freedom of the Drude oscillators to keep the atomic system close to the chosen temperature and the degrees of freedom of the DP at 1 K.36-37 The recent implementation of Thole atomic damping in the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)33 was used to reduce the short-range attraction of the respective Drude-core-atom interaction,


8

Page 9 of 31

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


where the parameter li is set to the suggested value of 2.6 Å for all simulations. The dielectric constant of 1,2-DCE is mainly determined by the variable molecular dipole moment and the resulting dipole moment distribution in the system. The reparametrization procedure for the development of the DO model thus needs to focus again on improving the molecular dipole moment by fitting the partial charges and improving the trans/gauche ratio by optimizing the torsional parameters.

2.2.1 Atomic polarizability A direct determination of the Drude charges via RESP fitting is not meaningful. Instead, the charge split is determined by means of the respective atomic polarizability α which is related to the force constant kD and the Drude charge qD by

α=

1 qD2 . 2 kD

The 1,2-DCE molecule is built from two equivalent CH2Cl subgroups, contributing equally the total polarizability. Not taking into account the hydrogen atoms as potential polarization sites, the DO model development in this study considered either the carbon or the chlorine of each CH2Cl subgroup as polarizable sites. It is worth noting that a DO model which assigned the polarizability to the chlorine atoms completely failed to reproduce the correct dielectric constant of bulk 1,2-DCE and exhibited a strong coupling between the DO degrees of freedom and the torsion. The experimental atomic polarizability α = 1.67 Å-3 was thus assigned to the carbon atom sites.41 The resulting partial charges for the different interaction sites of the 1,2-DCE-DO are collected in Table 1.

2.2.2 Torsion potential parameters A slightly modified approach to fit the dihedral parameter was applied to the polarizable model with the already fitted 1,2-DCE-NP model as a starting guess. It is important to mention that even though the potential energy profile for a single molecule might be reproduced correctly, the final conformational distribution can still vary to a large extent for the resulting condensed phase. Thus an additional


9


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

Page 10 of 31

modification of the dihedral parameters based on a set of MD simulations proved necessary for a proper reproduction of the conformational trans/gauche ratio. For the final set of parameters the obtained average dipole moment is close to the first-principles reference.11 The refitted dihedral parameters with the established partial charge set that can provide the proper reproduction of the conformational distribution in the bulk liquid simulation were selected and are listed in the Table 2.

2.3 Molecular dynamics simulation details All simulations were carried out by classical molecular dynamics using the LAMMPS simulation program package46. The intermolecular interaction potential used in this study is a Lennard-Jones (LJ) site-site potential plus Coulomb interactions. The 1/2 standard Lorentz-Berthelot rules, ε ij = (ε ii ε jj ) and σ ij = (σ ii + σ jj ) / 2 , were used

to derive the LJ potential parameters between unlike atom-types, where ε ij and σ ij are the energetic and size parameters of the LJ interaction between sites i and j, respectively. The nonbonded LJ interactions were treated using a regular spherical cut-off (cut-off radius rc = 10 Å) while the long-range Coulomb interaction was evaluated using the Ewald47-48 summation method with a precision of 10-6. During the entire simulation, periodic boundary conditions (PBC) were applied in all three directions and an integration time step of 1 fs was chosen. After initial energy minimization to remove unfavorable interactions, each system was equilibrated for 1 ns in the canonical (NVT) ensemble and subsequently for 3 ns in the isothermal-isobaric (NPT) ensemble at 1 atm and 298.15 K. Analysis programs were either taken from the Pizza.py toolkit49 or written in-house.

3. COMPUTATIONAL PROCEDURE After equilibration, production MD simulations were performed to compute thermodynamic and response properties for all optimized models of the CHCs CHCl3, CH2Cl2, 1,1-DCE, and 1,2-DCE. The properties evaluated were density, volumetric expansion coefficient, isothermal compressibility, heat capacity, dielectric constant,


10

Page 11 of 31

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


heat of vaporization, and surface tension. Two sets of simulations were thus required, corresponding to a “bulk liquid” system and a “liquid-vapor” system for each CHC. The remaining details for the bulk liquid and liquid-vapor simulation sets are given in the following paragraphs.

3.1 Simulation of bulk liquid systems For the equilibrium properties of each CHC, we prepared five independent bulk liquid systems. The lengths of the X, Y, Z edges of the cuboid simulation box are 22, 22, and 68 Å, respectively. The simulation cells contain 248 CHCl3, 312 CH2Cl2, 244 1,1-DCE, 252 1,2-DCE-NP, and 252 1,2-DCE-DO molecules to reproduce the bulk densities of the respective liquids at room temperature. The general settings for the bulk liquid simulations were given above.

3.1.1 Density ρ An average density ρ =

M was obtained from additional 1 ns constant pressure

simulations after the equilibration phase for each CHC. Densities were collected every 0.1 ps by monitoring the instantaneous volume V. Here and in the following, the denotes a statistical average over the different configurations of the system.

3.1.2 Volumetric expansion coefficient αP, isothermal compressibility kT, heat capacity CP and CV We used the fluctuations in the NPT ensemble to calculate the volumetric expansion coefficients αP, isothermal compressibility kT and heat capacity CP at constant pressure for all CHCs. For the calculations of αP, kT, and CP, additional NPT simulations were performed after equilibration at different temperatures T for all CHCs. Each simulation was run for 1 ns sampling the target properties every 0.1 ps during these production runs. For the computation of the target properties from fluctuations we followed Ref.50 using the following equations:


11


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

kT =

αP = CP =

Page 12 of 31

1 (< V 2 > − < V >2 ) NPT k BT < V > [< V ⋅U > − < V > ⋅ +P(< V 2 > − < V >2 )]NPT

,

k BT 2 < V > ik B N 1 + [( − 2 ) + 2P( − < V > ⋅ ) + P 2 (< V 2 > − < V >2 )]NPT 2 k BT 2

with temperature T, pressure P, volume V, Boltzmann’s constant kB, the potential energy U, the number of molecules N, and i, the number of degrees of freedom (DOF). Since none of the CHC molecules are restraint, we have i = 6.

To evaluate the constant volume heat capacity CV, additional NVT simulations were performed after equilibration at temperature T for all CHCs. Each system was simulated for 1 ns and samples were collected every 0.1 ps. The constant volume heat capacity CV was calculated from the fluctuations in the NVT ensemble: CV =

ik B N 1 + [( − 2 ) + 2P(< V ⋅U > − ⋅ ) + P 2 (< V 2 > − < V >2 )]NVT 2 k BT 2

3.1.3 Dielectric constant ε The static dielectric constant, ε, is related to the fluctuations of the system’s total dipole moment M in the computational box. Thus, the dielectric constant can be determined from the following fluctuation formula derived e.g. by Neumann51,

ε = ε∞ +

4π < M 2 > − < M > 2 . 3 Vk BT

Here, ε∞ denotes the high-frequency or optical dielectric constant34-35, which in the non-polarizable model is set to ε∞=1, V is the volume of the simulation box, kB is the Boltzmann constant and the system’s total dipole moment M (= Σµ) is obtained by summation over all molecular dipole moments in the simulation cell. As it has been observed before52, the fluctuation < M 2 > − < M > 2 tends to converge quite slowly (see below). For the NP force field simulations 12 ns long trajectories in the microcanonical (NVE) ensemble were required in addition to the equilibrations to ensure that the computed values are well converged. For the 1,2-DCE-DO simulation,


12

Page 13 of 31

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


the DO model in LAMMPS requires thermostats,33 as described above. Thus, a 12 ns simulation in the NVT ensemble was performed, treating the initial 2 ns as equilibration, with the following 10 ns used to record 10001 sample configurations at an interval of 1 ps for analysis.

3.2 Simulation of liquid-vapor interfacial systems The liquid-vapor interfacial systems are somewhat more complex to simulate in comparison to the bulk systems. Additionally, the surface tension is a property well-known to be rather problematic to reproduce with an NP force field.53-54 For the calculation of the liquid-vapor interfacial systems, five cubic simulation systems containing 1000 molecules each were prepared. The lengths of the cube edges of the respective simulation boxes are summarized in Table 3, which were selected to reproduce the same densities of the liquids as in the above bulk liquid simulations at room temperature. To obtain the liquid-vapor interfaces, the size of each simulation box in the Z direction was extended by a factor of 4, creating vacuum regions at the top and bottom side (PBC) of the CHC slabs thus generating a system with two liquid-vapor interfaces parallel to the XY-plane. Each system was simulated for 5 ns, with the initial 4 ns treated as the equilibration phase. During the final nanosecond samples were collected for the target properties every 0.1 ps. All simulations were run in the NVT ensemble using a larger cut-off radius rc = 15 Å. All other settings were consistent with the above simulation details.

Table 3. Sizes of the cubic liquid-vapor simulation boxes for the different CHC systems, only showing the edge length in the X dimension. Each system contained N=1000 molecules.

System

CHCl3

CH2Cl2

1,1-DCE

X (Å)

52.0

49.0

52.5

1,2-DCE-NP 1,2-DCE-DO


51.5

51.5

13


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

Page 14 of 31

3.2.1 Heat of vaporization △ Hvap The heat of vaporization △Hvap can be calculated from the difference between the potential energy of an ideal gas and the liquid state. Without considering additional corrections, the simplest approximation is given by ∆H vap (T ) ≈ (Hgas − Hliq) / N + RT , with Hgas and Hliq the potential energy of the N molecules in gas phase and liquid phase4, respectively, and R the gas constant. Hgas was taken from the potential energy of N molecules in a gas phase simulation while Hliq was deduced from the potential energy of each system in the liquid-vapor simulation.

3.2.2 Surface tension γ The surface tension was computed from the pressure tensor following the mechanical method of Kirkwood and Buff55 by integration of the difference between the normal pressure-tensor

component

Pn (z) = Pzz (z)

and

tangential

components

Pt (z) = (Pxx (z) + Pyy (z)) / 2 along the z direction. The resulting formula reads

γ (t) =

Lz  Px (t) + Py (t)   Pz (t) −  , where Pn is the pressure tensor in the direction n and 2 2 

Lz is the length of the simulation region along the z direction (perpendicular to the liquid-vapor interface).

4. RESULTS AND DISCUSSION In the following, the equilibrium properties, including density, volumetric expansion coefficient, isothermal compressibility, heat capacity, dielectric constant, heat of vaporization, and surface tension obtained from MD simulations performed as described above are presented. The models from the present work are labeled OPLS-NP and OPLS-DO to indicate either the standard nonpolarizable or polarizable Drude oscillator models reparametrized based on the original OPLS-AA force field. Table 4 compares the above physical properties of the CHCs computed with our new models with values obtained within the original OPLS-AA force field by Caleman et


14

Page 15 of 31

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


al.23,28, here labeled as OPLS, and also with experimental data40-42. As already mentioned, the statistics entering the dielectric constant are not easily converged and require simulations on the nanosecond time scale. We evaluated the dielectric constant as a function of simulation length for the CHC models (see next section). To control the influence of system size on the dielectric constant, we additionally calculated ε from three independent simulations with varying numbers of molecules for 1,2-DCE-NP and 1,2-DCE-DO, respectively. The results were identical to the values given in the table for the smaller systems. To extend the range of applicability of our force field model we also computed the isothermal compressibility and dielectric constant not only at room temperature but as a function of temperature for this class of the CHCs (see below).

Table 4. Properties of CHCs with the improved force field models. ρ, γ, ε, αP, κT, ΔHvap, CV, and CP denote the density, surface tension, relative dielectric constant, volumetric

expansion coefficient, isothermal compressibility, heat of vaporization, constant volume heat capacity and constant pressure heat capacity at 298.15 K. OPLS-NP is the new non-polarizable model which is directly comparable to the standard OPLS (all atom) parametrization. OPLS-DO denotes the polarizable model (see text for details). ρ Liquid CHCl3

CH2Cl2

1,1-DCE

∆Hvap

γ

αP

κT

CV

CP

(0.001/K)

(1/GPa)

(J/mol/K)

(J/mol/K)

ε

Model (g/cm3)

(kJ/mol)

(0.001N/m)

OPLS-NP

1.40

24.13

24.34

4.33

1.70

1.64

76.80

121.35

OPLS

1.37a

29.20a

11.8a

3.3a

2.28a

2.16a

89.6a

113a

Expt.

1.48d

31.28b

26.67c

4.71c

1.29b

1.03b

116.98b

117.00b

OPLS-NP

1.20

20.64

21.62

8.85

1.85

1.82

78.40

120.19

OPLS

1.21a

23.36a

10.5a

4.4a

2.31a

2.22a

75a

53.5a

Expt.

1.32b

28.82b

27.20c

8.82c

1.35b

1.03b

100.86b

100.88b

OPLS-NP

1.13

28.85

20.71

10.02

1.69

1.68

124.43

168.45

OPLS

1.18a

28.29a

15.1a

3.2a

1.60a

1.33a

100.7a

64a


15


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

1,2-DCE

a

1.17b

30.62b

24.07c

10.10c

1.33b

1.15b

126.18b

126.20b

OPLS-NP

1.16

29.66

27.57

10.79

1.31

1.15

126.84

164.96

OPLS-DO

1.21

31.31

29.95

14.01e

1.45

1.04

132.68

181.70

OPLS

1.24a

35.08a

23.2a

13.4a

1.43a

0.95a

106.7a

78a

Expt.

1.25b

35.16b

31.86c

10.13c

1.15b

0.82b

128.88b

128.90b

Ref.40

Ref.41

d e

Expt.

Ref.28

b c

Page 16 of 31

Ref.42

For the polarizable model ε∞=1.54 was determined according to the procedure from

Refs. 34-35 As Table 4 shows, the results obtained with our reparametrized models are in overall good agreement with the available experimental data. Our reparametrization leads either to an improved agreement with a range of experimental properties, or matches them roughly equally well as the previous OPLS-AA models. Additionally, it is important to stress that for the dielectric constant and surface tension, the reparametrized models reproduce the experimental value much more accurately than the original OPLS-AA force field. The improvement of the dielectric constants were achieved by assigning improved atomic partial charges yielding better agreement with the respective molecular dipole moments and, for the more flexible molecules, a more accurate conformational distribution obtained by reparametrization of the torsion potential coefficients. We obtained 32.3% and 30.8% of trans conformer for 1,2-DCE-NP and 1,2-DCE-DO model, respectively, which is consistent with the experimental value of 34.9%56 and in good agreement with an ab-initio molecular dynamics calculation (32.2%)11. Both new models clearly outperform the original OPLS-AA force field in this respect, which results in 9.7%. The most notable properties which are not improved over the standard OPLS-AA and are still off of the experimental values in the NP models are the density and the heat of vaporization, which might be the price to pay for improving the dielectric performance. There


16

Page 17 of 31

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


seems to be a clear trade-off for what can be achieved with a NP model, where we opted for a rather balanced treatment of important solvent properties. Some improvement w.r.t. the underestimated density might also be possible using larger simulations thus reducing statistical fluctuations and possible PBC effects. It is also noted that nonideal gas corrections in the estimates of the heat of vaporization might improve the agreement with experiment. In both cases we adopted a more conservative approach (e.g. traditional approximation in 3.2.1 for △Hvap) to render our results more easily comparable with former parametrizations.

The introduction of polarizable sites, thus increasing the flexibility of the respective model, could improve the density and thermal compressibility since both long-range attraction and short-range repulsion are influenced by the additional Coulomb potential terms.57 The density ρ and thermal compressibility kT obtained within the 1,2-DCE-DO model are indeed noticeably improved over the 1,2-DCE-NP model, illustrating the expected correlation of these two properties and confirming the assumption that the introduction of polarizability would yield improved ρ and kT. For the liquid-vapor interfacial simulations of 1,2-DCE, the Drude oscillator model also leads to a significant improvement of agreement with experiments for quantities △Hvap and γ, compared to 1,2-DCE-NP model without loosing the agreement for other properties. The heat capacities CP and Cv are hardly influenced, showing that no energy is artificially deposited in the degrees of freedom modelling the electronic polarizability. The dielectric constant ε of 1,2-DCE obtained with the OPLS-DO model is in general agreement with the experimental reference, however not as good as with the OPLS-NP model. Since the system is assumed to be isotropic, < M > 2 ≈ 0 , the overestimated value of ε according to equation in 3.1.3 implies that the term

< M 2 > /V of the polarizable system is increased. This can be explained by the fact that the dipole moment fluctuations in the numerator of (< M 2 > − < M > 2 ) / V increase owing to the broadened gauche dipole moment regime while at the same


17


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

Page 18 of 31

time also the density increases.

4.1 Dielectric constant as a function of simulation time and system size We evaluated the effects of different simulation parameters such as the trajectory length and system size during the determination of the static dielectric constant at 298.15 K. Figure 1 shows the estimates for the dielectric constant as a function of sampling time for the different CHCs.

Figure 1. Static bulk dielectric constant ε as a function of sampling time for the CHCs studied. All the simulations had already been equilibrated prior to these runs.

From the Figure 1, it is clear that a properly converged value can only be achieved after about 6 ns. In general, the fluctuation properties are more difficult to predict than simple averages, which means that longer equilibration and production simulation durations are needed. For the current work we settled on simulations with 10 ns production time length to obtain consistent estimates of dielectric constant.


18

Page 19 of 31

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


To study the dependence of the computed dielectric constant on the system size, we focused on the 1,2-DCE models, since this system exhibits the most subtle effects due to its internal degree of freedom which dominates the molecular dipole moment. As Figure 2 shows, no monotonous dependence of the dielectric constants for both nonpolarizable and polarizable 1,2-DCE molecules on system size can be established. This result agrees well with Gereben’s work52, but is unfortunately different from Morrow’s finding58 that the static dielectric constant increases with increasing system size and slowly converges to the theoretical of infinite system size within PBC. Preliminary results in a different context seem to indicate that much larger system sizes might be required to enter the regime of smoother convergence and the systematic coarse-graining of models on a molecular scale to yield dielectric continuum models definitely needs to be studied more thoroughly.52 For the purpose of this reparametrization we settled on the larger systems with the number of molecules stated above in section 3.1.

Figure 2. Static bulk dielectric constant ε as a function of system size for 1,2-DCE-NP and 1,2-DCE-DO models. The number of molecules N1, N2, N3, N4 in the four


19


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

Page 20 of 31

different simulation systems are 63, 126, 189 and 252, respectively.

4.2 Temperature dependence of isothermal compressibility (kT) and dielectric constant (ε) We evaluated the isothermal compressibility and static dielectric constant as a function of temperature for the CHCs to investigate how well our reparametrized models capture the temperature dependence in the equilibrium bulk over a temperature range of relevance in applications. Additional simulations were run at temperatures from 283.15 to 328.15 K, in 5 K increments (see Figures 3 and 4), where temperatures above 320 K are already beyond the boiling point of the bulk liquid CH2Cl2 and thus not included for this molecule. All other settings were consistent with the above bulk liquid simulation details. The isothermal compressibility and dielectric constant are fitted to a 2nd-order polynomial in T temperatures (as it is also used in the Handbook of Chemistry and Physics41) over the above range of temperatures in order to permit interpolation and regularization. The resulting coefficients are given in Table 5 and 6 for isothermal compressibility and dielectric constant, respectively. The interpolation polynomials were also used in order to compare the simulations to experimental data. Corresponding plots of the isothermal compressibility kT(T) and dielectric constant ε(T) as a function of temperature and their fits for each CHC with respect to the experimental results are presented in Figure S1 and S2 in the Supporting Information.


20

Page 21 of 31

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


Figure 3. Temperature dependence of kT for the CHC models. The black dashed lines are the polynomial fitting curves for each liquid shown. Outliers, which have been removed from the fit, are indicated by a cyan circle. These are for the CHCl3 model (red dot) at 328 K, for the CH2Cl2 model (green triangle) at 303 K and for the 1,2-DCE-NP model (black cross) at 303 K and 318 K. All models capture the trend of increasing kT with increasing temperature.

Table 5. Parametrization of the temperature dependence of isothermal compressibility constants in polynomial form kT = A + BT + CT2. N is the number of data points included in each fit (see Figure 3). Molecules

N

Tmin

Tmax

A

B

C

CHCl3

9

283.15

328.15

9.319e+00

-7.238e-02

1.548e-04

CH2Cl2

7

283.15

318.15*

-2.949e+00

5.964e-03

3.305e-05

1,1-DCE

10

283.15

328.15

9.164e+00

-7.035e-02

1.517e-04

1,2-DCE-NP

8

283.15

328.15

7.727e+00

-5.679e-02

1.173e-04

1,2-DCE-DO

10

283.15

328.15

1.228e+00

-9.403e-03

2.948e-05


21


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

Page 22 of 31

* Tmax for CH2Cl2 is limited by the boiling point.

Figure 4. Temperature dependence of the static bulk dielectric constant ε for the CHC models. The black dashed lines are the polynomial fitting curves for each liquid shown. All models capture the trend of decreasing ε with increasing temperature.

Table 6. Parametrization of temperature dependence of dielectric constants in polynomial form ε = A + BT + CT2. N is the number of data points included in each fit. Molecules

N

Tmin

Tmax

A

B

C

CHCl3

10

273.15

328.15

5.627e+00

1.551e-02

-6.606e-05

CH2Cl2

8

273.15

318.15*

4.717e+01

-1.975e-01

2.292e-04

1,1-DCE

10

273.15

328.15

3.061e+01

-8.209e-02

4.260e-05

1,2-DCE-NP

10

273.15

328.15

6.764e+01

-3.437e-01

5.146e-04

1,2-DCE-DO

10

273.15

328.15

2.582e+02 -1.529e+00

2.396e-03

* Tmax for CH2Cl2 is limited by the boiling point.


22

Page 23 of 31

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


As Figures 3 and 4 show, within the range of relevant temperatures, all models exhibit the same physically reasonable behavior: kT(T) increases with increasing temperature while ε(T) decreases with the increase of temperature, which is also in agreement with the experimental values (see Support Information) and the original OPLS-AA force field.28 It has to be noted that all models not only capture the trend of ε(T) but, with the exception of CH2Cl2 close to its boiling point, also reproduce the experimental slope of this trend rather well over a range of 40K according to the plots in Figure S1 in the Supporting Information.

5. CONCLUSIONS In the present work, MD simulations were performed to assess the validity of reparametrized OPLS-AA based force field models in predicting solvent properties of a set of important CHCs. The new parameter sets yield a range of physical properties of these CHCs in excellent agreement with the available experimental data at ambient conditions. Additionally, it is possible to reproduce the temperature dependence of the isothermal compressibility and dielectric constant. The nonpolarizable force field reparametrizations are based on the original OPLS-AA force field as starting point. The atomic partial charges and thereby the molecular dipole moments were improved via the RESP38-39 approach based on first-principle computations. Systematically improved torsional potential energy surface cuts for the internal degree of freedom in the flexible 1,1- and 1,2-DCE molecules were constructed by adjusting the torsional potential parameters to match corresponding first-principle calculations and cross-validated based on rotational conformer distributions from short MD simulations. This systematic approach which combines carefully adjusted atomic charges, yielding a better model for the respective molecular dipole moments, with a more accurate representation of the conformational distribution, obtained by reparametrization of the torsion potential coefficients, was able to improve significantly on earlier models for the dielectric constant without sacrificing other solvent properties. For the 1,2-DCE molecule, we also developed a polarizable Drude oscillator (DO) model to better capture the more intricate response properties of this


23


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

Page 24 of 31

molecule whose molecular dipole moment is determined by the torsional degree of freedom. The most promising model for 1,2-DCE assigns polarizability to the carbon atom site. A careful modification of the torsion potential coefficients starting from the nonpolarizable parametrization proved to be crucial, which allows the DO model to remedy some problems with correlations among physical properties observed in the reparametrization of the NP force field.

This work presents a simple but effective systematically first-principles based re-fitting protocol for the derivation of improved parameter sets for CHC MD models to reproduce accurate dielectric constants without sacrificing well-reproduced physical properties of existing models. The methodology developed here might be useful in deriving and improving further force fields for (organic) solvents relevant in applications that are sensitive to response properties.

ASSOCIATED CONTENT Supporting Information The Supporting information is available free of charge on the ACS publication website at DOI: Table S1 and S2 in the supporting information show the values of the isothermal compressibility and dielectric constant for each CHC at different temperatures. Temperature dependence plots of the isothermal compressibility and dielectric constant and their fits for each CHC with respect to the experimental results are presented in Figure S1 and S2.

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected] Notes The authors declare no competing financial interest.


24

Page 25 of 31

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


ACKNOWLEDGEMENTS We gratefully acknowledge support by the International Centre for Energy Research (ICER), a cooperation project between Technical University of Munich (TUM) and Nanyang Technological University (NTU) in Singapore. Computational resources have been partially provided by HPC projects at the Leibniz Rechenzentrum (LRZ), Garching. Z. L. is thankful for the support from the China Scholarship Council.

REFERENCES

(1) Wood, J. M., Chlorinated Hydrocarbons: Oxidation in the Biosphere. Environ. Sci. Technol. 1982, 16, 291A-297A. (2) Dang, L. X., Intermolecular Interactions of Liquid Dichloromethane and Equilibrium Properties of Liquid–Vapor and Liquid–Liquid Interfaces: A Molecular Dynamics Study. J. Chem. Phys. 1999, 110, 10113-10122. (3) Lopes, P. E.; Harder, E.; Roux, B.; Mackerell Jr, A. D., Formalisms for the Explicit Inclusion of Electronic Polarizability in Molecular Modeling and Dynamics Studies. In Multi-Scale Quantum Models for Biocatalysis, Springer: 2009; pp 219-257. (4) Fennell, C. J.; Li, L.; Dill, K. A., Simple Liquid Models with Corrected Dielectric Constants. J. Phys. Chem. B 2012, 116, 6936-6944. (5) Benjamin, I., Theoretical Study of the Water/1,2 ‐ Dichloroethane Interface: Structure, Dynamics, and Conformational Equilibria at the Liquid–Liquid Interface. J. Chem. Phys. 1992, 97, 1432-1445. (6) Millot, C.; Rivail, J., Molecular Dynamics Simulation of Conformational Equilibrium of 1,2-Dichloroethane. J. Mol. Liq. 1989, 43, 1-11. (7) Salas, F. J.; Méndez-Maldonado, G. A.; Núñez-Rojas, E.; Aguilar-Pineda, G. E.; Dominguez, H.; Alejandre, J., Systematic Procedure to Parametrize Force Fields for Molecular Fluids. J. Chem. Theory Comput. 2015, 11, 683-693.


25


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

Page 26 of 31

(8) Jorgensen, W. L.; Binning Jr, R. C.; Bigot, B., Structures and Properties of Organic Liquids: N-Butane and 1,2-Dichloroethane and Their Conformation Equilibriums. J. Am. Chem. Soc. 1981, 103, 4393-4399. (9) Taniguchi, Y.; Takaya, H.; Wong, P.; Whalley, E., Effect of Pressure on Molecular Conformations. Ii. Trans–Gauche Equilibrium of 1,2 ‐ Dichloroethane and 1,2 ‐ Dibromoethane. J. Chem. Phys. 1981, 75, 4815-4822. (10) Kuo, K.-H.; Lin, Y.-S.; Shih, J.-J.; Liao, Y.-H.; Lin, J.-S.; Yang, C.-M.; Lin, J.-L., 1,2-Dichloroethane on Cu (100) and Oxygen-Covered Cu (100): Conformation, Orientation, and Transformation. J. Phys. Chem. C 2007, 111, 7757-7764. (11) Murugan, N. A.; Hugosson, H. W.; Ågren, H., Solvent Dependence on Conformational

Transition,

Dipole

Moment,

and

Molecular

Geometry

of

1,2-Dichloroethane: Insight from Car− Parrinello Molecular Dynamics Calculations. J. Phys. Chem. B 2008, 112, 14673-14677. (12) Schlaich, A.; Knapp, E. W.; Netz, R. R., Water Dielectric Effects in Planar Confinement. Phys. Rev. Lett. 2016, 117, 048001. (13) Bonthuis, D. J.; Netz, R. R., Beyond the Continuum: How Molecular Solvent Structure Affects Electrostatics and Hydrodynamics at Solid–Electrolyte Interfaces. J. Phys. Chem. B 2013, 117, 11397-11413. (14) Bonthuis, D. J.; Gekle, S.; Netz, R. R., Profile of the Static Permittivity Tensor of Water at Interfaces: Consequences for Capacitance, Hydration Interaction and Ion Adsorption. Langmuir 2012, 28, 7679-7694. (15) Bonthuis, D. J.; Gekle, S.; Netz, R. R., Dielectric Profile of Interfacial Water and Its Effect on Double-Layer Capacitance. Phys. Rev. Lett. 2011, 107, 166102. (16) Ringe, S.; Oberhofer, H.; Reuter, K., Transferable Ionic Parameters for First-Principles Poisson-Boltzmann Solvation Calculations: Neutral Solutes in Aqueous Monovalent Salt Solutions. J. Chem. Phys. 2017, 146, 134103. (17) Ringe,

S.;

Oberhofer,

H.;

Hille,

C.;

Matera,

S.;

Reuter,

K.,

Function-Space-Based Solution Scheme for the Size-Modified Poisson–Boltzmann Equation in Full-Potential Dft. J. Chem. Theory Comput. 2016, 12, 4052-4066.


26

Page 27 of 31

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


(18) Berger, D.; Logsdail, A. J.; Oberhofer, H.; Farrow, M. R.; Catlow, C. R. A.; Sherwood, P.; Sokol, A. A.; Blum, V.; Reuter, K., Embedded-Cluster Calculations in a Numeric Atomic Orbital Density-Functional Theory Framework. J. Chem. Phys. 2014, 141, 024105. (19) Sánchez, V. M.; Sued, M.; Scherlis, D. A., First-Principles Molecular Dynamics Simulations at Solid-Liquid Interfaces with a Continuum Solvent. J. Chem. Phys. 2009, 131, 174108. (20) Iozzi, M. F.; Cossi, M.; Improta, R.; Rega, N.; Barone, V., A Polarizable Continuum Approach for the Study of Heterogeneous Dielectric Environments. J. Chem. Phys. 2006, 124, 184103. (21) Jorgensen, W. L.; Tirado-Rives, J., Potential Energy Functions for Atomic-Level Simulations of Water and Organic and Biomolecular Systems. Proc. Nati. Acad. Sci. U. S. A. 2005, 102, 6665-6670. (22) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A., Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157-1174. (23) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D., Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant. J. Chem. Theory Comput. 2011, 8, 61-74. (24) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I., Charmm General Force Field: A Force Field for Drug ‐ Like Molecules Compatible with the Charmm All ‐ Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671-690. (25) Jorgensen, W. L.; Ibrahim, M., Structures and Properties of Organic Liquids: N-Alkyl Ethers and Their Conformational Equilibriums. J. Am. Chem. Soc. 1981, 103, 3976-3985. (26) Jorgensen, W. L.; McDonald, N. A.; Selmi, M.; Rablen, P. R., Importance of


27


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

Page 28 of 31

Polarization for Dipolar Solutes in Low-Dielectric Media: 1,2-Dichloroethane and Water in Cyclohexane. J. Am. Chem. Soc. 1995, 117, 11809-11810. (27) 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. (28) van der Spoel, D.; van Maaren, P. J.; Caleman, C., Gromacs Molecule & Liquid Database. Bioinformatics 2012, 28, 752-753. (29) Fox, T.; Kollman, P. A., Application of the RESP Methodology in the Parametrization of Organic Solvents. J. Phys. Chem. B 1998, 102, 8070-8079. (30) Chen, S.; Yi, S.; Gao, W.; Zuo, C.; Hu, Z., Force Field Development for Organic Molecules: Modifying Dihedral and 1–N Pair Interaction Parameters. J. Comput. Chem. 2015, 36, 376-384. (31) Dang, L. X.; Chang, T.-M., Molecular Dynamics Study of Water Clusters, Liquid, and Liquid–Vapor Interface of Water with Many-Body Potentials. J. Chem. Phys. 1997, 106, 8149-8159. (32) Lamoureux, G.; Roux, B., Modeling Induced Polarization with Classical Drude Oscillators: Theory and Molecular Dynamics Simulation Algorithm. J. Chem. Phys. 2003, 119, 3025-3039. (33) Dequidt, A.; Devémy, J.; Pádua, A. A., Thermalized Drude Oscillators with the Lammps Molecular Dynamics Simulator. J. Chem. Inf. Model. 2015, 56, 260-268. (34) Lamoureux, G.; MacKerell Jr, A. D.; Roux, B., A Simple Polarizable Model of Water Based on Classical Drude Oscillators. J. Chem. Phys. 2003, 119, 5185-5197. (35) Vorobyov, I. V.; Anisimov, V. M.; MacKerell, A. D., Polarizable Empirical Force Field for Alkanes Based on the Classical Drude Oscillator Model. J. Phys. Chem. B 2005, 109, 18988-18999. (36) Sprik, M.; Klein, M. L., A Polarizable Model for Water Using Distributed Charge Sites. J. Chem. Phys. 1988, 89, 7556-7560. (37) Van Belle, D.; Wodak, S., Extended Lagrangian Formalism Applied to Temperature Control and Electronic Polarization Effects in Molecular Dynamics Simulations. Comput. Phys. Commun. 1995, 91, 253-262.


28

Page 29 of 31

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


(38) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A., A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269-10280. (39) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollmann, P. A., Application of RESP Charges to Calculate Conformational Energies, Hydrogen Bond Energies, and Free Energies of Solvation. J. Am. Chem. Soc. 1993, 115, 9620-9631. (40) Marcus, Y., The Properties of Solvents: Wiley Series in Solutions Chemistry. John Wiley & Sons Ltd: West Sussex, UK: 1998. (41) Haynes, W. M., CRC Handbook of Chemistry and Physics; CRC press, 2014. (42) Frenkel, M.; Hong, X.; Dong, Q.; Yan, X.; Chirico, R., 2.2. 1 C1-C2. In Densities of Halohydrocarbons, Springer: 2003; pp 107-130. (43) Frisch, M.; Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Montgomery, J.; Vreven, T.; Kudin, K.; Burant, J., Gaussian 03, Revision C. 02. 2008. (44) Inspyred: Bio-inspired Algorithms in Python — inspyred 1.0 Documentation. http://inspyred.github.io/. Accessed: 6 May 2014. (45) Dunning Jr, T. H., Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007-1023. (46) Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1-19. (47) 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. (48) 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. (49) Plimpton, S.; Jons, M., Pizza.py Toolkit. http://pizza.sandia.gov/. 2005. (50) Allen, M. P.; Tildesley, D. J., Computer Simulation of Liquids; Oxford university press, 1989. (51) Neumann, M., Dipole Moment Fluctuation Formulas in Computer Simulations of Polar Systems. Mol. Phys. 1983, 50, 841-858. (52) Gereben, O.; Pusztai, L., On the Accurate Calculation of the Dielectric Constant


29


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

Page 30 of 31

from Molecular Dynamics Simulations: The Case of SPC/E and SWM4-DP Water. Chem. Phys. Lett. 2011, 507, 80-83. (53) Caleman, C.; Hub, J. S.; van Maaren, P. J.; van der Spoel, D., Atomistic Simulation of Ion Solvation in Water Explains Surface Preference of Halides. Proc. Nati. Acad. Sci. U. S. A. 2011, 108, 6838-6842. (54) Wensink, E. J.; Hoffmann, A. C.; van Maaren, P. J.; van der Spoel, D., Dynamic Properties of Water/Alcohol Mixtures Studied by Computer Simulation. J. Chem. Phys. 2003, 119, 7308-7317. (55) Kirkwood, J. G.; Buff, F. P., The Statistical Mechanical Theory of Surface Tension. J. Chem. Phys. 1949, 17, 338-343. (56) Tanabe, K., Calculation of Infrared Band Intensities and Determination of Energy Differences of Rotational Isomers of 1,2-Dichloro-, 1,2-Dibromo- and 1-Chloro-2-Bromoethane. Spectrochim. Acta Mol. Biomol. Spectrosc. 1972, 28, 407-424. (57) Bedrov, D.; Borodin, O.; Li, Z.; Smith, G. D., Influence of Polarization on Structural, Thermodynamic, and Dynamic Properties of Ionic Liquids Obtained from Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 4984-4997. (58) Morrow, T.; Smith, E., Simulation Calculation of Dielectric Constants: Comparison of Methods on an Exactly Solvable Model. J. Stat. Phys. 1990, 61, 187-201.


30

Page 31 of 31

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


TOC Graphic


31

Benchmarks and Dielectric Constants for Reparametrized OPLS and

Recommend Documents