Applicability of Effective Fragment Potential Version 2-Molecular

Dec 10, 2018 - Effective fragment potential version 2-molecular dynamics (EFP2-MD) simulations, where the EFP2 is a polarizable force field based on a...
1 downloads 0 Views 771KB Size
Subscriber access provided by AUSTRALIAN NATIONAL UNIV

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

Applicability of Effective Fragment Potential Version 2 Molecular Dynamics (EFP2-MD) Simulations for Predicting Dynamic Liquid Properties Including Supercritical Fluid Phase Nahoko Kuroki, and Hirotoshi Mori J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b07446 • Publication Date (Web): 10 Dec 2018 Downloaded from http://pubs.acs.org on December 16, 2018

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 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 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.

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 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

The Journal of Physical Chemistry

Applicability of Effective Fragment Potential Version 2 - Molecular Dynamics (EFP2-MD) Simulations for Predicting Dynamic Liquid Properties Including Supercritical Fluid Phase Nahoko Kuroki1 and Hirotoshi Mori2,* 1Department

of Chemistry and Biochemistry, Graduate School of Humanities and Sciences,

Ochanomizu University, 2-1-1 Otsuka, Bunkyo-ku, Tokyo 112-8610, Japan 2Faculty

of Core Research Natural Science Division, Ochanomizu University, 2-1-1 Otsuka, Bunkyo-ku, Tokyo 112-8610, Japan

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 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 19

ABSTRACT Effective fragment potential version 2 - molecular dynamics (EFP2-MD) simulations, where the EFP2 is a polarizable force field based on ab initio electronic structure calculations, were applied to predict the static and dynamic liquid properties of compressed liquid NH3. By analyzing the temperature dependency of the radial distribution function, the auto-correlation functions of velocity (Cv(t)) and reorientation (Cr(t)), and the self-diffusion constant, we clarified that the ab initio EFP2 force field can effectively describe the properties of compressed liquids. These descriptions can be performed with at least semi-quantitative accuracy and at a sufficiently low computational cost. In the EFP2-MD protocol, no force field training is required. This training is mandatory when simulating liquid properties with classical MD techniques (especially in extreme conditions with high pressures and temperatures). EFP2MD is a promising technique for predicting the physicochemical properties of novel functional compressed liquids, including super critical fluid phase properties.

ACS Paragon Plus Environment

2

Page 3 of 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

The Journal of Physical Chemistry

1. Introduction Molecular-level insights into the transporting properties of liquid materials (particularly under extreme conditions involving high pressures and temperatures) are important in the fields of chemical engineering and physical solution chemistry.1 For example, supercritical fluids made up of small organic molecules are used both as extraction solvents and as novel organic synthesis environments.2,3 Supercritical fluids can be mixed with the other supercritical fluids or gases because intermolecular interactions among their constituent molecules are weak, i.e. their diffusion constants are very large. With the backgrounds, predicting the physicochemical nature of functional compressed liquids (including supercritical fluids) has attracted the interest of many engineers. The self-diffusion constant is one of the most important dynamic properties of liquid materials. Researchers have used classical molecular dynamics (MD) simulations to evaluate the self-diffusion constants of functional liquids under various extreme conditions.4-10 As is well known, however, classical MD simulations are based on classical force fields, which describe the intra/intermolecular interaction potential of the target systems. Generally, classical force fields are optimized for reproducing 1) experimentally observed target material properties or 2) intermolecular interactions evaluated by high-level quantum chemical calculations for various molecular configurations. Without great effort, it is very difficult to obtain sufficiently accurate classical force field parameters for novel (unknown) functional liquid materials especially under extreme conditions.

In order to give predictive ability of the functional liquid

materials to the MD simulations more precisely, it is immediately urged to apply first principle MD (AIMD) calculations (e.g. Car-Parrinello MD; CPMD11,12, fragment molecular orbital based MD; FMOMD13-18, etc.), in which intermolecular interaction is evaluated in quantum mechanical level. Because the computational costs of AIMD are much higher than the computational costs of classical MDs, it has been difficult in computational chemistry to perform AIMD simulations to predict novel functions of compressed liquid materials, including supercritical fluids. Recently, we showed that it is possible to predict mixed solvent properties using effective fragment potential version 2 - MD simulations (EFP2-MD).19 In the EFP2-MD protocol, internal ACS Paragon Plus Environment

3

The Journal of Physical Chemistry 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 19

degrees of freedom in the target molecules are fixed, we focus only on relative molecular motions. The intermolecular potential is evaluated by describing each interaction term (i.e., electrostatic, exchangerepulsion, polarization, and dispersion interaction) by compactly pre-expanding each monomer wave function (For details, see the next section).20,21 Because the EFP2-MD method is based on a quantum mechanical wave function for monomer fragments, it is possible to apply the method to various liquid systems as a relatively cheap AIMD simulation with predictive accuracy.19,22-25 EFP2-MD has wellbalanced characteristics of low computational cost and semi-quantitative chemical accuracy. These characteristics make EFP2-MD promising for applications involving the prediction of the physicochemical nature of novel compressed liquids. In this paper, we discuss the applicability of EFP2-MD simulations to the prediction of static and dynamic liquid properties under high pressure conditions, including the properties of supercritical fluids. As a test system, we focus on highly compressed liquid NH3.

2. Details of Effective fragment potential version 2 (EFP2) theory Here, we explain EFP2 theory in more detail. The original version of the EFP, which is often called as EFP1, was developed by Day and Gordon et al.20,26 Focusing in only hydrating systems, in the EFP1 theory, intermolecular interaction in the target systems were tried to describe with compactly expanded localized wave function of constituent monomer fragment molecules, in which repulsive interaction term was modeled with some empirical parameters. The EFP2 is an upgrade theory of the original EFP1 to suit for not only water solvent but also general solvents without using any empirical parameters.20,21 In the EFP2 theory, total inter-fragment interaction is evaluated as the sum of the four energy components. 𝐸TOTAL = 𝐸ES + 𝐸EXREP + 𝐸POL + 𝐸DISP

(1)

, where EES, EEXREP, EPOL and EDISP are electrostatic, exchange-repulsion, polarization, and dispersion energy, respectively. ACS Paragon Plus Environment

4

Page 5 of 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

The Journal of Physical Chemistry

Electrostatic interaction is calculated by multipole expansion, of which expansion points are at the atoms and bond-midpoints generated by means of Stone distributed multipolar analysis (DMA).27 In the case of short range, Coulomb interactions are damped in consideration of the fragments overlap to avoid excessive instability.28-30 Exchange-repulsion is evaluated as an intermolecular overlap truncated at the quadratic term.31-33 Kinetic and overlap integrals are calculated between each pair of fragments on the fly using wave function for each fragment.

EXREP

𝐸

=

∑ 𝑘,𝑙

{

―2ln |𝑆𝑘𝑙|𝑆𝑘𝑙2 ―4 ― 2𝑆𝑘𝑙 𝐹A𝑘𝑖𝑆𝑖𝑙 + 𝐹B𝑙𝑗𝑆𝑗𝑘 ― 2𝑇𝑘𝑙 𝑅 𝜋 𝑘𝑙 𝑖∈A 𝑗∈B 𝑍𝐽 𝑍𝐼 1 1 1 +2𝑆𝑘𝑙2 ― ― +2 +2 ― 𝑅𝑘𝐽 𝑅𝑙𝐼 𝑅𝑘𝑗 𝑅𝑖𝑙 𝑅𝑘𝑙

(∑

(

}

)











𝐽∈B

𝐼∈A

𝑗∈B

𝑖∈A

)

(2) , where i, j, k, and l are localized molecular orbitals; I and J are nuclei, S and T are the overlap and kinetic energy integrals, respectively; and F is the intramolecular Fock matrix. Polarization interaction is treated as induced dipoles of one fragment with static multipole field and the induced field of the other fragment. The truncated polarizability term is evaluated as a tensor sum of anisotropic localized molecular orbitals polarizabilities.34 Dispersion interaction, which is an interaction among the induced dipoles, is expressed by an inverse R expansion. 𝐸DISP =

𝐶𝑛

∑𝑅

𝑛

𝑛 = 6,8

(3) , where the 1/R6 term is so called van der Waals interaction term, and the 1/R8 term is added in order to correct short range interaction.35 As seen in the above equations, EFP2 theory is based on well-defined molecular fragmentation with no chemical reaction. In EFP2-MD simulation, the EFP2 were applied to describe intermolecular potentials in an ab initio manner.

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry 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 19

3. Computational Details The purpose of this study is to assess whether EFP2-MD simulations can predict compressed liquids properties.

First, the molecular structure of NH3 was optimized using quantum chemical

calculations at the MP2/aug-cc-pVTZ level of theory.36 With respect to the optimum structures, an ab initio polarizable EFP2 force field was determined using the “MAKEFP” module implemented in the GAMESS-US (Aug. 2016 ver.) quantum chemical program package.37

Figure 1

Coordinate system of NH3 dimer.

Before performing MD simulations, we evaluated EFP2 force field accuracy by comparing the intermolecular interaction potential surfaces obtained by a set of EFP2 and quantum chemical calculations at the CCSD(T)/aug-cc-pVTZ level of theory. The structural parameters used in the potential surface evaluation are shown in Figure 1. To compare EFP2 and CCSD(T) (in terms of both their total interaction energies and their physicochemical components, i.e. electrostatic, exchange repulsion, polarization, and dispersion components), the CCSD(T) interaction energies were decomposed using a localized molecular orbital energy decomposition analysis (LMO-EDA) method.38 In the LMO-EDA calculations, we applied a counterpoise method to correct basis set superposition errors.

ACS Paragon Plus Environment

6

Page 7 of 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

The Journal of Physical Chemistry

After the EFP2 accuracy evaluation, we subsequently performed a set of EFP2-MD simulations for the super critical and compressed liquid phases of NH3. In the EFP2-MD simulations, we used a set of cubic periodic boxes containing 500 NH3 molecules having an isothermal-isobaric (NPT) ensemble, with a cutoff distance of 10 Å. The time step was set to 1 fs/step for the EFP2-MD simulations. The pressure condition was set to a constant value of 500 bar. While the temperature conditions were varied as 203, 243, 323, 373, 423, 473, and 700 K, we theoretically evaluated the density and diffusion constant of NH3. The Nosé-Hoover and Hoover’s equations were applied for thermal and pressure controls, respectively. With the above conditions, a set of 0.8 ns equilibration and 1.0 ns production runs were performed to evaluate the radial distribution functions (RDFs), autocorrelation velocity functions (Cv(t)), reorientation functions (Cr(t)), and self-diffusion constants.

4. Results and Discussion 4-1. Accuracy of the EFP2 force field Figure 2 depicts intermolecular interaction potential energy surfaces for NH3 dimers (at R = 3.5 and φ = 0; see definitions of R and φ in Figure 2). The figure also gives the mean error (in kcal mol-1) for the EFP2 force field from the highly accurate quantum chemical CCSD(T)/aug-cc-pVTZ level of theory. Figure 2 clearly shows that the EFP2 force field follows the corresponding CCSD(T) results within chemical accuracy. If we examine the potential data in more detail, EFP2 overestimates instability around the structure of (θ1,θ2)=(150,40) relative to the CCSD(T) results. This occurs because of the EFP2’s tendency to underestimate dispersion interaction, which originates from the fact that the EFP2 dispersion interaction is expanded by the Hartree-Fock wave function of interacting monomer fragments. However, the difference in interaction potential between EFP2 and CCSD(T) is less than 1.5 kcal mol-1. This means that the EFP2 force field has sufficient accuracy for evaluating intermolecular interactions in the condensed phase of NH3.

ACS Paragon Plus Environment

7

stable

Page 8 of 19

-0.5 ± 0.4

in cm-1

-0.1 ± 0.1

TOTAL

0.1 ± 0.1

DISP

0.1 ± 0.2

POL+CT

Figure 2

EFP2

Mean Error in kcal mol-1

-0.4 ± 0.1

EXREP ES

ab initio

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

unstable

The Journal of Physical Chemistry

Two-dimensional potential energy surfaces obtained by LMO-EDA and EFP2 methods. For the LMO-EDA, we applied an ab initio CCSD(T)/aug-cc-pVTZ level of quantum chemical calculation.

ACS Paragon Plus Environment

8

Page 9 of 19

4-2. Static and dynamic properties of compressed/supercritical NH3 described with EFP2-MD simulations A set of N-N RDFs for different temperature conditions are depicted in Figure 3. To obtain the RDFs, we statistically averaged trajectory data along 150 ps EFP2-MD simulations. Comparing with the previous classical MD report using specially parametrized force fields for reproducing intermolecular interaction potential by quantum chemistry calculations7, our RDFs gave basically similar structures. In the RDFs, the NH3 distributions (N atom) around other N atoms are given for the temperature conditions. Under relatively low temperature conditions, we observed a set of apparent peaks assigned to the first and second NH3 hydration shells. Increasing temperatures, the second RDF peak vanishes, resulting in a very broad band. With the temperature ascension, we also observed a decrease in the NH3 coordination number (nNN), which is evaluated as the integral of the RDF along the radial coordinate (See Table 1). These observations indicate that increasing temperatures lead to a phase transition from a compressed liquid to a supercritical fluid.

2.5

203K 243K 323K 373K 423K 473K 700K

2.0

gNN(r)

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

1.5 1.0 0.5 0.0

0

2

4

6

8

10

r (Å)

Figure 3

Temperature dependence of N-N radial distribution functions (RDFs) for compressed NH3.

ACS Paragon Plus Environment

9

The Journal of Physical Chemistry 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

Page 10 of 19

Temperature dependencies of NH3 coordination numbers (nNN(r)) for the first solvation shell (in the region of r = 2.0–5.0 Å). T/K

nNN(r)

203

11.6

243

10.6

323

8.2

373

6.6

423

5.2

473

4.2

700

2.2

The liquid-supercritical NH3 densities in the temperature range of 203-700K were 0.741-0.164 g cm-3 and 0.652-0.134 g cm-3 by experimental7 and EFP2-MD, respectively. EFP2-MD gave ca. -10% systematic error. It should be reminded that rigid rotor approximation was applied in the EFP2 theory. Thus, theoretical densities obtained by the EFP2-MD should be smaller than the experimental values, since all molecules in the system have frozen degrees of freedom. On the other hand, the corresponding density change by OPLS3-MD was 0.820-0.143 g cm-3.

That is, OPLS3-MD gave ca. +10%

overestimation and -10% underestimation with varying temperature condition. In the case of density prediction for novel liquids including supercritical fluids, EFP2-MD with a set of systematic error that can be explained by proper physicochemical reasons would be preferred. NH3 confined in a high-pressure box is polarized, because the average intermolecular distance is shorter than it is under physicochemically mild conditions. Figure 4 depicts variations in the temporal dipole moment of NH3 simulated by EFP2-MD simulations (T = 203, 473 and 700 K). We observed that almost all NH3 molecules were strongly polarized in the compressed liquid phase (T = 203 K). The average variation in dipole moment was +0.65 Debye under these conditions, meaning that the strong polarization of NH3 must be considered. When heating the system to induce a transition into the supercritical fluid phase (T = 473 and 700 K), the distribution of the strongly polarized NH3 was observed to be sparse with the elongation of averaged intermolecular distance. It should be noted that ACS Paragon Plus Environment

10

Page 11 of 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

The Journal of Physical Chemistry

several NH3 molecules in supercritical fluid phase remained strongly polarized. This can be explained by the fact that microscopic density fluctuates in supercritical fluids (i.e., several NH3 molecules form micro clusters). When higher temperature conditions were simulated, the population of strongly polarized NH3 molecules decreased.

Figure 4

Temporal dipole moment variations for NH3 simulated by EFP2-MD simulations (T = 203, 473 and 700 K). The variations were defined as differences from the dipole moment of free NH3 molecules predicted by EFP2 theory (1.61 Debye, which agrees well with experimental data of 1.47 Debye.39) For clarity, NH3 molecules are depicted as single atoms placed at their center of mass, and a set of representative snapshot structures with time evolution are depicted.

The temperature dependencies of the auto-correlation functions for velocity (Cv(t)) and reorientation (Cr(t)) are given in Figure 5. Defining u as a unit vector along the principal C3 axis of NH3, Cr(t) can be represented by the following equation: 𝐶𝑟(𝑡) = 〈𝒖𝑖(𝑡)𝒖𝑖(0)〉 = 〈𝑐𝑜𝑠𝜃(𝑡)〉

(4),

where (t) is the angle formed by the principal C3 axes of NH3 at t = 0 and t. From Figure 5(a), we observed that the rate of Cr(t) decay increases with increasing temperatures. This can be well explained by the fact that liquid density decreases with increasing temperatures. In low density conditions, Cr(t) diminishes more rapidly, because molecular rotation is easy. The temperature dependency of Cr(t) simulated by our EFP2-MD method appropriately followed previous classical MD reports by Vyalov.40 The following equation defines the auto-correlation function for velocity (Cv(t)):

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry 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

〈𝒗𝑖(𝑡)𝒗𝑖(0)〉

𝐶𝑣(𝑡) = 〈𝒗 (0)2〉 𝑖

Page 12 of 19

(5)

Cv(t) is a function that describes the relative translational motion of the target molecule in their condensed phase. For an ideal gas, Cv(t) is always equal to 1, because of the lack of molecular collision. As shown in Figure 5(b), a set of exponential decays of Cv(t) is observed for relatively high temperature (i.e. low density) conditions (323-700 K). In contrast, according to our simulations, Cv(t) becomes negative for t = 0.1–0.5 ps in the compressed liquid phase at low temperature (203 K). The negative value of Cv(t) means that the translation direction is scattered by molecular collision, because of the high density of the compressed liquid. Analyzing dynamic molecular motion, it was clarified that we successfully modeled two different condensed phases (the compressed liquid phase and the supercritical fluid phase), using an EFP2-MD simulation with changing temperature conditions.

Figure 5

Temperature dependence of Cr(t) and Cv(t) for compressed NH3.

Finally, it would be useful to demonstrate the merits of applying the EFP2-MD scheme to further study liquid materials. In Figure 6, the self-diffusion constants of NH3 predicted by the EFP2MD simulation were shown as functions of temperature and liquid density. In this study, the selfdiffusion constants (D) were defined according to the following equation, using the long-time limits of the mean-square displacement (MSD).

ACS Paragon Plus Environment

12

Page 13 of 19

The Journal of Physical Chemistry 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

𝐷 = lim 6𝑡〈|𝒓𝑖(𝑡) ― 𝒓𝑖(0)|2〉

(6)

𝑡→∞

As shown in Fig. 6, classical MD with optimized FF, which is specially optimized for reproducing intermolecular interaction potential by quantum chemistry calculations, well followed the experimental self-diffusion constants by Groß et al.41 On the other hand, EFP2-MD and OPLS3-MD42 gave slight overestimates of self-diffusion constants compared with the experimental results. It should be noted that our object in this study is not to reproduce experimental data in an exact manner but to predict dynamic physicochemical property such as self-diffusion constant prior to performing experiments with semi quantitative accuracy and low computational costs. Generally speaking, evaluation of diffusion constants is difficult even to reproduce their order without applying well trained force fields. Fig. 6 shows that EFP2-MD has succeeded in reproducing the order of diffusion constant with no effort for optimizing such an empirical parameter. The overall results presented in this report indicate that a new stage of functional compressed liquid design (including supercritical fluid phase) can be developed with EFP2-MD, a simple and sufficiently accurate ab initio MD scheme. It is well known that a set of long term MD simulations (at least sub-ns order) is mandatory for predicting dynamic molecular properties like self-diffusion constants. Although there is significant need in the chemical engineering field to perform ab initio MD simulations for predicting liquid functionality (particularly for compressed liquids including supercritical fluids), it is currently very difficult to perform such predictive simulations, because of their high computational costs relative to classical simulations.

Thus, classical MD

simulations with well-trained force field parameters have been used for engineering applications. Considering the balance between accuracy and computational and simulation setup costs, EFP2-MD is a promising protocol for designing functional liquids (including the physicochemical properties of such liquids) under high pressures and temperatures. The computational cost of EFP2-MD to track 1 ns trajectory for the target system is half week with modern in-house computers (32 CPU cores), which is much faster than so called ab initio MD (e.g. CPMD, FMO-MD, etc.). In our group, physicochemical property predictions for various supercritical fluids are now in progress. ACS Paragon Plus Environment

13

The Journal of Physical Chemistry 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 6

Page 14 of 19

Dependency of self-diffusion constant of NH3 on temperature. Experimental and classical MD (optimized FF) data were those reported by Orabi et al.7

5. Concluding Remarks In this paper, the applicability of simulating strongly compressed liquids using EFP2-MD, for including supercritical fluid phase, was reported for NH3 using a set of numerical examples. While the EFP2 force field was constructed based only on the ab initio wave function for a monomer (NH3) fragment, we showed that both static (RDFs) and dynamic condensed phase properties (auto correlation functions and self-diffusion constants) could be evaluated with sufficient accuracy compared with the previous experimental and theoretical data. No empirical parameterization, which is mandatory in classical MD simulations, was required before performing MD simulations with the EFP2-MD protocol, even for compressed liquids and supercritical fluids. Overall, we showed that it is possible to predict the dynamic liquid properties of compressed liquids at least semi-quantitatively using the simple EFP2MD protocol described in this report.

AUTHOR INFORMATION ACS Paragon Plus Environment

14

Page 15 of 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

The Journal of Physical Chemistry

Corresponding Author *E-mail: [email protected]. Tel.: +81-3-5978-5718. Fax: +81-3-5978-5717. ORCID Nahoko Kuroki: 0000-0002-2615-3481 Hirotoshi Mori: 0000-0002-7662-9636 Notes The authors declare no competing financial interest.

ACKNOWLEDGMENTS Some of the presented calculations were performed at the Research Center for Computational Science (RCCS), the Okazaki Research Facilities, and the National Institutes of Natural Sciences (NINS).

This study was supported in part by the Advanced Information and Communication

Technology for Innovation (ACT-I; Grant number: JPMJPR16UB) Precursory Research for Embryonic Science and Technology programs (PRESTO; Grant number: JPMJPR16NC) from the Japan Science and Technology (JST) Agency, and JSPS KAKENHI (Grant Number: 16K13928 and 18J11490).

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry 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 16 of 19

REFERENCES (1) Hendriks, E.; Kontogeorgis, G. M.; Dohrn, R.; de Hemptinne, J. -C.; Economou, I. G.; Zilnik, L. F.; Vesovic, V. Industrial Requirements for Thermodynamics and Transport Properties. Ind. Eng. Chem. Res. 2010, 49(22), 11131-11141. (2) Vyhmeister, E.; Valdes-Gonzalez, H.; Muscat, A. J.; Suleiman, D.; Estevez, L. A. Surface Modification of Porous Silicon-Based Films Using Dichlorosilanes Dissolved in Supercritical Carbon Dioxide. Ind. Eng. Chem. Res. 2013, 52(13), 4762-4771. (3) Tomita, K.; Machmudah, S.; Quitain, A. T.; Sasaki, M.; Fukuzato, R.; Goto, M. Extraction and Solubility Evaluation of Functional Seed Oil in Supercritical Carbon Dioxide. J. Supercrit. Fluids 2013, 79, 109-113. (4) Feng, H.; Liu, X.; Gao, W.; Chen, X.; Wang, J.; Chen, L.; Lüdemann, H. -D. Evolution of SelfDiffusion and Local Structure in Some Amines over a Wide Temperature Range at High Pressures: A Molecular Dynamics Simulation Study. Phys. Chem. Chem. Phys. 2010, 12(45), 15007-15017. (5) Idrissi, A.; Vyalov, I.; Kiselev, M.; Fedorov, M. V.; Jedlovszky, P. Heterogeneity of the Local Structure in Sub- and Supercritical Ammonia: A Voronoi Polyhedra Analysis. J. Phys. Chem. B 2011, 115(31), 9646-9652. (6) Guevara-Carrion, G.; Vrabec, J.; Hasse, H. Prediction of Transport Properties of Liquid Ammonia and Its Binary Mixture with Methanol by Molecular Simulation. Int. J. Thermophys. 2012, 33(3), 449-468. (7) Orabi, E. A.; Lamoureux, G. Polarizable Interaction Model for Liquid, Supercritical, and Aqueous Ammonia. J. Chem. Theory Comput. 2013, 9(4), 2035-2051. (8) Alberti, M.; Amat, A.; Farrera, L.; Pirani, F. From the (NH3)2-5 Clusters to Liquid Ammonia: Molecular Dynamics Simulations Using the NVE and NpT Ensembles. J. Mol. Liq. 2015, 212, 307315. (9) Ushiki, I.; Ota, M.; Sato, Y.; Inomata, H. A Kinetic Study of Organic Compounds (Acetone, Toluene, n-Hexane and n-Decane) Adsorption Behavior on Activated Carbon Under Supercritical Carbon Dioxide Conditions at Temperature from 313 to 353 K and at Pressure from 4.2 to 15.0 MPa. J. Supercrit. Fluids 2014, 95, 187-194. (10) Sohrevardi, N.; Bozorgmehr, M. R.; Heravi, M. M.; Khanpour, M. Transport Properties of Mixtures Composed of Acetone, Water, and Supercritical Carbon Dioxide by Molecular Dynamics Simulation. J. Supercrit. Fluids 2017, 130, 321-326. (11) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471-2474. (12) Saharay, M.; Balasubramanian, S. Electron Donor-Acceptor Interactions in Ethanol-CO2 Mixtures: An Ab Initio Molecular Dynamics Study of Supercritical Carbon Dioxide. J. Phys. Chem. B 2006, 110(8), 3782-3790. (13) Kato, Y.; Fujiwara, T.; Komeiji, Y.; Nakano, T.; Mori, H.; Okiyama, Y.; Mochizuki, Y. Fragment Molecular Orbital-Based Molecular Dynamics (FMO-MD) Simulations on Hydrated Cu(II) Ion. Chem-Bio Inf. J. 2014, 14, 1-13.

ACS Paragon Plus Environment

16

Page 17 of 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

The Journal of Physical Chemistry

(14) Matsuda, A.; Mori, H. Theoretical Study on the Hydration Structure of Divalent Radium Ion Using Fragment Molecular Orbital-Molecular Dynamics (FMO-MD) Simulation. J. Solution Chem. 2014, 43(9), 1669-1675. (15) Mori, H.; Hirayama, N.; Komeiji, Y.; Mochizuki, Y. Differences in Hydration Between cis- and trans-Platin: Quantum Insights by Ab Initio Fragment Molecular Orbital-Based Molecular Dynamics (FMO-MD). Comput. Theor. Chem. 2012, 986, 30-34. (16) Komeiji, Y.; Mochizuki, Y.; Nakano, T.; Mori, H. Recent Advances in Fragment Molecular Orbital-Based Molecular Dynamics (FMO-MD) Simulations in Molecular Dynamics; InTech: Rijeka, Croatia, 2012. (17) Fujiwara, T.; Mori, H.; Mochizuki, Y.; Osanai, Y.; Miyoshi, E. 4f-in-Core Model Core Potentials for Trivalent Lanthanides. Chem. Phys. Lett. 2011, 510(4), 261-266. (18) Fujiwara, T.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Mori, H.; Nakano, T.; Miyoshi, E. Fragment Molecular Orbital-Based Molecular Dynamics (FMO-MD) Simulations on Hydrated Zn(II) Ion. Chem. Phys. Lett. 2010, 490(1), 41-45. (19) Kuroki, N.; Mori, H. Applicability of Effective Fragment Potential Version 2 - Molecular Dynamics (EFP2-MD) Simulations for Predicting Excess Properties of Mixed Solvents. Chem. Phys. Lett. 2018, 694, 82-85. (20) Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An Effective Fragment Method for Modeling Solvent Effects in Quantum Mechanical Calculations. J. Chem. Phys. 1996, 105, 1968-1986. (21) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112(1), 632-672. (22) Kuroki, N.; Mori, H. Effective Fragment Potential Version 2 - Molecular Dynamics (EFP2-MD) Simulation for Investigating Solution Structures of Ionic Liquids. Chem. Lett. 2016, 45(8), 10091011. (23) Hands, M. D.; Slipchenko, L. V. Intermolecular Interactions in Complex Liquids: Effective Fragment Potential Investigation of Water-tert-Butanol Mixtures. J. Phys. Chem. B 2012, 116(9), 2775-2786. (24) Ghosh, M. K.; Cho, S. G.; Choi, C. H. A Priori Prediction of Heats of Vaporization and Sublimation by EFP2-MD. J. Phys. Chem. B 2014, 118(18), 4876-4882. (25) Ghosh, M. K.; Cho, S. G.; Choi, T. H.; Choi, C. H. A Priori Predictions of Molecular Density by EFP2-MD. Theor. Chem. Acc. 2016, 135(12), 254/1-254/7. (26) Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. The Effective Fragment Potential Method: A QM-Based MM Approach to Modeling Environmental Effects in Chemistry. J. Phys. Chem. A 2001, 105(2), 293-307. (27) Stone, A. J. Distributed Multipole Analysis, or How to Describe a Molecular Charge Distribution. Chem. Phys. Lett. 1981, 83(2), 233-239. (28) Slipchenko, L. V.; Gordon, M. S. Electrostatic Energy in the Effective Fragment Potential Method: Theory and Application to Benzene Dimer. J. Comput. Chem. 2007, 28(1), 276-291.

ACS Paragon Plus Environment

17

The Journal of Physical Chemistry 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 19

(29) Freitag, M. A.; Gordon, M. S.; Jensen, J. H.; Stevens, W. J. Evaluation of Charge Penetration Between Distributed Multipolar Expansions. J. Chem. Phys. 2000, 112(17), 7300-7306. (30) Slipchenko, L. V.; Gordon, M. S. Damping Functions in the Effective Fragment Potential Method. Mol. Phys. 2009, 107(8-12), 999-1016. (31) Jensen, J. H.; Gordon, M. S. An Approximate Formula for the Intermolecular Pauli Repulsion Between Closed Shell Molecules. Mol. Phys. 1996, 89(5), 1313-1325. (32) Jensen, J. H.; Gordon, M. S. An Approximate Formula for the Intermolecular Pauli Repulsion Between Closed Shell Molecules. II. Application to the Effective Fragment Potential Method. J. Chem. Phys. 1998, 108(12), 4772-4782. (33) Jensen, J. H. Modeling Intermolecular Exchange Integrals Between Nonorthogonal Molecular Orbitals. J. Chem. Phys. 1996, 104(19), 7795-7796. (34) Li, H.; Netzloff, H. M.; Gordon, M. S. Gradients of the Polarization Energy in the Effective Fragment Potential Method. J. Chem. Phys. 2006, 125(19), 194103/1-194103/9. (35) Adamovic, I.; Gordon, M. S. Dynamic Polarizability, Dispersion Coefficient C6 and Dispersion Energy in the Effective Fragment Potential Method. Mol. Phys. 2005, 103(2-3), 379-387. (36) 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(2), 1007-1023. (37) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S., et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14(11), 1347-1363. (38) Su, P.; Li, H. Energy Decomposition Analysis of Covalent Bonds and Intermolecular Interactions. J. Chem. Phys. 2009, 131(1), 014102/1-014102/15. (39) Nelson, Jr. R. D.; Lide, D. R.; Maryott, A. A. Selected Values of electric dipole moments for molecules in the gas phase; NBS10; NSRDS: Washington, D. C., USA, 1967. (40) Vyalov, I.; Kiselev, M.; Tassaing, T.; Soetens, J. C.; Federov, M.; Damay, P.; Idrissi, A. Reorientation Relaxation in Supercritical Ammonia. J. Mol. Liq. 2011, 159(1), 31-37. (41) Groß, T.; Buchhauser, J.; Price, W. E.; Tarassov, I. N.; Lüdemann, H. -D. The p,T-Dependence of Self-Diffusion in Fluid Ammonia. J. Mol. Liq. 1997, 73-74, 433-444. (42) Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight, J. L., et al. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins. J. Chem. Theory Comput. 2016, 12(1), 281-296.

ACS Paragon Plus Environment

18

Page 19 of 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

The Journal of Physical Chemistry

TOC graphic

ACS Paragon Plus Environment

19