Subscriber access provided by University of South Dakota
C: Surfaces, Interfaces, Porous Materials, and Catalysis
Development of Reaction Density Functional Theory and Its Application to Glycine Tautomerization Reaction in Aqueous Solution Weiqiang Tang, Cheng Cai, Shuangliang Zhao, and Honglai Liu J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b05383 • Publication Date (Web): 22 Aug 2018 Downloaded from http://pubs.acs.org on August 23, 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 35 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
Development of Reaction Density Functional Theory and Its Application to
2
Glycine Tautomerization Reaction in Aqueous Solution
3
Weiqiang Tang1, Cheng Cai2, Shuangliang Zhao1,1 and Honglai Liu2 1
4
State Key Laboratory of Chemical Engineering and School of Chemical Engineering, East
5
China University of Science and Technology, Shanghai, 200237, China 2
6
School of Chemistry and Molecular Engineering, East China University of Science and
7
Technology, Shanghai, 200237, China
8
Abstract
9
Accurate determination of reaction barrier is crucial in chemical reaction engineering since
10
it’s directly associated with the reaction rate and reaction path. Whereas massive quantum
11
calculations relying on commercial software have been reported for understanding reaction
12
mechanisms, the solvent effect on the reaction barrier is still poorly understood. The difficulty
13
originates partially from the ill capture of inhomogeneous solvent structure around the
14
reaction site. Herein, we propose a reaction density functional theory (RxDFT) by combining
15
quantum density functional theory for calculating intrinsic reaction energy with classical
16
density functional theory for addressing solvation contribution. For demonstration, the
17
glycine tautomerization reaction in aqueous solution is revisited by using RxDFT. We
18
illustrate that the free energy profiles of two distinct reaction paths can be well predicted, in
19
comparison with the results of MD/MP2 simulation. Since no solvent molecule is explicitly
20
considered, the developed RxDFT is promising to provide an efficient alternative to QM/MM
21
approach or ab initio calculation.
22
1
To whom correspondence should be addressed. Email:
[email protected] 1
ACS Paragon Plus Environment
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
1
I.
Introduction
2
It has been well recognized that reaction solvent has a profound effect on the reaction rate,
3
chemical equilibrium and even the reaction mechanism.1-3 The use of solvents in organic
4
chemical reactions, especially in synthesis, has become a regulation method. Using different
5
solvents, the reaction efficiencies are quite different.4 For example, in the substitution reaction
6
of butyl bromide and sodium cyanide: C 4 H 9 Br + NaCN → C 4 H 9 CN + NaBr , the solvent plays
7
a pivotal role.5 Namely, if butyl bromide is simply mixed with sodium cyanide in aqueous
8
solution, the reaction does not occur even if refluxed at 373 K for two weeks. This is mainly
9
due to the fact that butyl bromide is insoluble in water and the substrate cannot contact with the
10
reagents. In contrast, the reaction can take place in the presence of alcoholic solvents although
11
the reaction rate and yield are still quite low. Interestingly, the reaction rate in
12
dimethylformamide (DMF) solution is ten times faster than that in alcohol. The influences of
13
the solvent on the reaction process are collectively referred to as the solvent effect or the
14
medium effect of the reaction.4
15
Generally, there are three categories of methods for addressing the solvent effects of
16
chemical reactions. The first one is the zeroth-order approximation of solvent effect,2-3 i.e.,
17
solvent is treated as continuum, and the solvent effect is reflected by the parameter setting of
18
dielectric constant during the calculation. In other words, the reaction is considered to occur in
19
a vacuum or gas phase or in a continuum. At present, most of the quantum chemistry
20
calculations are based on this approach.6-7 This method has the advantage of low computational
21
cost and can make reasonable explanation for gaseous phase reactions. However, in most
22
liquid phase reactions, especially when the solvent is composed of polar or charged molecules
23
and in addition the reactant and product molecules have polarity changes during the reaction 2
ACS Paragon Plus Environment
Page 2 of 35
Page 3 of 35 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
process, there is a significant and even qualitative difference in the calculated reaction free
2
energy with or without properly considering the solvent’s role.1, 8 This discrepancy is certainly
3
crucial to determine the reaction path and reaction mechanism.
4
The second category refers to the art of simulations coupling classical mechanics and
5
quantum mechanics (QM/MM), firstly proposed by Car and Parrinello, which takes the solvent
6
effect strictly into account through multiscale calculations.9 In principle, the reaction kinetics
7
including the reaction path can be exactly described by using QM/MM. The essential idea of
8
QM/MM method is to divide the reaction system into two regions.9-10 The core zone (quantum
9
system) undergoes an intrinsic reaction, which generally contains only the molecules and ions
10
involved in the reaction (sometimes including catalysts). The peripheral zone is a solvent
11
system that did not participate in the reaction, but interacts with the quantum system to affect
12
the occurrence of the reaction and the choice of the reaction path. Reaction substrates that have
13
not reacted yet are also included in the solvent system. The solvent system can be described by
14
molecular mechanics, whose influence on the quantum system is reflected in the crossover
15
Hamiltonian operator.10
16
In general, QM/MM calculation is theoretically complicated and computationally
17
expensive,11 and thus its extensive application has been limited. The zeroth-order
18
approximation may bring qualitative errors. Therefore, researchers introduce several more
19
advanced statistical mechanics theories to describe the solvent effect by combined quantum
20
mechanics and statistical mechanics.12-16 So far, a series of approaches have been developed,
21
including multi-configuration self-consistent field theory (MCSCF),17-18 the averaged solvent
22
electrostatic potential theory (ASEP),19-21 mean field approximation (MFA),22 and the 3
ACS Paragon Plus Environment
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
1
reference interacting site model self-consistent field (RISM/SCF) method.23 The most notable
2
method probably is the RISM/SCF theory which was developed based on RISM. RISM
3
represents an intermediate approach between the continuum description and molecular
4
simulation method, in which no solvent molecules are explicitly considered while the solvent
5
structure is captured by means of the solvent density distribution function.24-27 The solvent
6
density distribution can be obtained by solving the RISM equation iteratively coupled with a
7
closure relation. Within the framework of RISM/SCF, the solvent structure is affected by the
8
chemical nature and molecular configuration of solute (reaction subsystem) through the
9
solute-solvent interaction. Once the solute (namely, the chemical nature and/or molecular
10
configuration) changes during a reaction, the solvent structure is updated by solving the RISM
11
equation again. The system reaches the final self-consistent equilibrium between the reaction
12
sub-system and the solvent after multiple iterations.28
13
Whereas the RISM/SCF method can describe the solvent effect in a self-consistent
14
manner with reasonable computational cost,29 its extension to investigating the reaction
15
dynamics is restricted. This is because the involved RISM theory belongs to the scope of
16
integral equation methods, which are basically established through the correlation functions
17
of equilibrium systems. In other words, RISM is not tractable for describing dynamic
18
properties such as diffusion.
19
Similar to RISM theory, classical density functional theory (CDFT) represents an
20
efficient and robust modern statistical mechanics theory.30 Unlike quantum density functional
21
theory (QDFT), which relies on electron density distribution as functional variable, CDFT
22
takes the spatial density distribution of molecular species as functional variable, and it allows 4
ACS Paragon Plus Environment
Page 4 of 35
Page 5 of 35 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
us to predict the structure and thermodynamic properties of large-scale fluid systems involving
2
no chemical reaction.31 Depending on the geometrical complexity of system components,
3
three different versions of CDFT have been developed: atomic DFT, molecular DFT and
4
polymeric DFT.30 As reflected by their names, atomic DFT can be employed to investigate
5
simple fluid systems composed of spherical particles,22 while molecular DFT can be adopted
6
to investigate molecular fluid systems by taking into account of the orientational contribution
7
of non-spherical rigid molecules.32-33 In a further step, the chain-connectivity contribution in
8
polymeric molecules is considered in polymeric DFT.34-35 Although the functionals in these
9
three versions of CDFT are notably different, they share the same variational framework.30
10
Besides these sub-branches of CDFT, a special variant of molecular DFT, i.e., site DFT,
11
should be mentioned as well. Site DFT, originally proposed by Chandler et al.,36 can be
12
utilized to study molecular fluid systems by dealing with the molecules with their atomic site
13
densities.33, 37 However, limited by the constraint of rigid-molecular geometry, the site density
14
in site DFT is not an independent variable.33, 38
15
In recent decades, CDFT has been applied to tackle various problems ranging from
16
fundamental physics to engineering challenges. To name but a few, CDFT was extended for
17
accurate prediction of solvation free energy, and recently it has been adopted for fast
18
screening of hydration free energy.39-40 In parallel, the screening of gas adsorption in porous
19
materials was also conducted for targeting porous materials with high gas storage
20
capacity.41-42 The extension of CDFT to non-equilibrium field is also successful. The
21
dynamical version of atomic DFT, conventionally called dynamic DFT, has been extensively
22
developed and widely applied in the computational-aided design of supercapacitor with high
23
performance.43-46
5
ACS Paragon Plus Environment
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
1
In light of the success and robustness of CDFT, an interesting question arises on the
2
combination of CDFT with QDFT for assessing the solvent effects of chemical reactions. The
3
first attempt was made probably by Arias and his co-workers,47 who proposed a joint DFT
4
very recently. Within the framework of joint DFT for a solute-solvent system, the solvent is
5
commonly treated with the polarizable continuum model, in which the solvent density is
6
directly related with the solute electron density, rather than an independent variable.38, 48 In
7
addition, the communication between solute and solvent is described with the so-called
8
density-only electronic den7sity functional method,48-49 and this makes the electron density to
9
be the unique variable in the joint DFT, allowing for robust and efficient minimization
10
procedure.50 Although explicit description of water system was individually developed later
11
on along the scope of site DFT by the same group, the deviations of the predicted oxygen and
12
hydrogen density profiles from the experimental results are apparent,38 and furthermore the
13
incorporation of this extended site DFT into the joint DFT has been rarely reported.
14
Herein we report the first attempt on the combination of QDFT and molecular DFT in
15
different length scales to investigate the chemical reaction in aqueous solution. By using this
16
multiscale description, we propose a reaction density functional theory (RxDFT). To
17
demonstrate, the RxDFT is thereafter applied to describe the geometry and thermodynamics of
18
the tautomerization reaction of glycine in aqueous solution. The remainder of this work is
19
organized as follows: In next section, the theoretical framework of RxDFT is introduced,
20
along with the illustrative computational details for the tautomerization reaction. The results
21
are discussed in the sec. III, followed by a short conclusion in sec. IV.
22
II.
23
1.
Modelling and Theory Theoretical framework of RxDFT
6
ACS Paragon Plus Environment
Page 6 of 35
Page 7 of 35 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 2
Figure 1. Thermodynamic circle for a representative chemical reaction (A+BC) occurs in
3
gas phase and in solution. FsolA , FsolB and FsolC are the solvation free energies of reagent and
4
product molecules, respectively, and E R denotes the reaction free energy in gas phase, and
5
FR is the reaction free energy in solution.
6 7
Considering a representative simple chemical reaction in a solvent: A+BC, where A
8
and B denote two different reagent molecules, and C denotes the product molecule. To
9
address the solvent effect, a thermodynamic circle can be established between the reaction in
10
dilute gas phase and that in the solvent, as depicted in Fig. 1. Since the solvation free energy
11
gives the reversible minimum work of moving a solute molecule from dilute gas phase to the
12
solvent,39 the solvation process of molecules A, B or C in the chemical reaction can be
13
C characterized with the individual solvation free energies, namely, FsolA , FsolB and Fsol . The
14
reaction free energy is denoted as FR in the solvent and E R in the gas phase. Note that in
15
dilute gas phase, the reaction free energy is usually approximated with the reaction energy as
16
the entropic contribution is overall neglectable. 7
ACS Paragon Plus Environment
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
1 2 3 4 5
Page 8 of 35
Following the thermodynamic circle, the reaction free energy in the solvent can be calculated by:
FR = ER + FsolC − FsolA − FsolB .
(1)
By setting ∆FS = FsolC − FsolA − FsolB , the above equation can be further simplified as:
FR = ER + ∆FS .
(2)
6
Here ∆FS represents the change in the solvation free energy after the occurrence of
7
chemical reaction. Clearly, when ∆FS is relatively small, eq.(2) recovers to the conventional
8
zeroth-order approximation of the solvent effect. However, if ∆FS is not ignorable, the
9
contribution of reaction free energy originates from two parts, as similar to the RISM/SCF
10
method,1 namely, the intrinsic reaction free energy and the solvation free energy.
11
Correspondingly, we may divide the system into two regions: the core region (quantum
12
system) undergoing an intrinsic reaction, and the peripheral region (solvent system). Whereas
13
several interesting attempts have been reported regarding the choice of the interface between
14
the core region and peripheral region,51-54 e.g., the solvent molecules in the first layer
15
surrounding the intrinsic reaction can be included into the core region; here we take the core
16
region as the one containing only the molecules and ions involved in the reaction.
17
The solvation process influences the mechanisms of chemical reactions in solution, but
18
the relaxation of the solvent is often much faster than the reactive event.4, 55-56 Namely,
19
solvation free energy calculation can be decoupled with the computation of the intrinsic
20
reaction free energy. While the intrinsic reaction free energy can be individually computed
21
with QDFT, the solvation free energy can be addressed with CDFT, being expressed as a
22
functional of the solvent density distribution function surrounding the solute: 8
ACS Paragon Plus Environment
Page 9 of 35 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
Fsol [ ρ ( x);Vext (x); µ , V , T ] = Θ[ ρ ( x);Vext (x); µ , V , T ] − Θ[ ρ 0 ;Vext (x); µ , V , T ] ,
2
where ρ (x) is the density profile of solvent molecules with x=(r, Ω) refers to the coordinates
3
and spatial orientation of solvent molecules. ρ 0 and µ are the solvent bulk density and bulk
4
chemical potential, respectively. V and T are the system volume and absolute temperature,
5
and they are given in prior and thus omitted below in the set of functional variables. Vext ( x )
6
represents the external potential to the solvent system, and it originates from the
7
solvent-solute interaction. Clearly, Vext ( x) changes when the molecular geometry and
8
charge distribution of the quantum system change duo to the occurrence of chemical reaction.
9
Certainly, Vext ( x ) is responsible for the interaction between the quantum system and solvent
10
system, and it comprises typically the coulomb interaction and Lennard-Jones interaction, as
11
disposed in other multiscale approaches.23, 47, 57 At a given Vext ( x ) , the grand potential can be
12
calculated as:
(3)
13
Θ[ ρ ( x);Vext ( x)] = k BT ∫ dx ρ ( x ){ln[ ρ ( x) Λ 3 ] − 1} + F exc [ ρ ( x )] + ∫ dx[Vext ( x) − µ ]ρ ( x ) ,
14
where k B is the Boltzmann constant; Λ is an effective thermal wavelength of solvent
15
molecule. F ex [ ρ ( x )] is the excess free energy functional accounting for the solvent
16
intermolecular interaction, of which no exact expression is available. To date, several
17
treatments have been developed for formulating the excess free energy functional. The
18
detailed formulation of this term will be discussed below.
(4)
19
The solvent density distribution minimizes the free energy functional of the solvent
20
system under the thermodynamic equilibrium condition. In other words, the solvent density
9
ACS Paragon Plus Environment
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 10 of 35
1
distribution and the solvation free energy of the solute can be obtained by minimizing
2
solvation free energy functional in eq.(3),39, 58 namely:
3
δ Fsol [ ρ (x);Vext ( x)] =0 ⇒ δρ (x)
ρ ( x ) = ρ eq ( x) .
(5)
F = Fsol
4
Because CDFT directly relates the free energy to the density distribution of system
5
components and bypasses the redundant thermodynamic integral, it greatly reduces the
6
computational cost and displays great superiority compared with molecular simulations.58
7
With the help of CDFT, the hydration free energies of different neutral solutes and ions can
8
be quantitatively predicted in an efficient way.39 On the other hand, compared with the RISM
9
theory, the theoretical framework of CDFT is more flexible, and it has recently been extended
10 11
to non-equilibrium system to describe the diffusion dynamics. 2.
Reaction free energy and activation free energy
12
The evaluation of activation free energy and reaction free energy represents one of the
13
most important routes to explore reaction mechanisms. From the perspective of physical
14
chemistry, the reaction free energy supplies an important criterion for judging the pathway of
15
reaction, and the magnitude of the activation energy determines the reaction probability.
16
According to the transition state theory,59 as shown in Fig. 2a, the free energies of the quantum
17
system in gas phase before and after reaction are denoted as F0 and F2 , respectively.
18
Conventionally, the reaction free energy is calculated as ER = F2 − F0 . Assuming that the free
19
energy of transition state system is F1 , the energy barrier or activation free energy of reaction
20
is E A = F1 − F0 .
10
ACS Paragon Plus Environment
Page 11 of 35 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 2
Figure 2. Typical energy profiles for a reaction a) in gas phase, and b) in a solvent. The
3
difference between the energies of transition state and reactant is the energy barrier or
4
activation free energy of reaction; it represents the minimum energy that a reacting system
5
must acquire for the transformation to take place. The difference between the energies of the
6
product and the reactant is the free energy of reaction.
7 8
Similarly, the reaction free energy is calculated as FR = F2′ − F0 in a solution, as shown
9
in the Fig. 2b. The energy barrier or activation free energy of reaction is FA = F1′ − F0 . Owning
10
to the decoupling between the intrinsic reaction free energy and solvent effect, FR
11
FA )
12
reaction with the solvation free energy difference, as stated in Eq.(2).
13
3. Revisit of the tautomerization reaction of glycine
(and
can be evaluated by summing up the corresponding E R (and E A ) for the intrinsic
14
Tautomerization is pervasive in chemistry and biology.60 There have been growing
15
theoretical interests in studying tautomerization in solution,57, 61-69 and this is because not
16
only
tautomerization reactions in solution can provide tractable model systems,70-72 but also
11
ACS Paragon Plus Environment
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 12 of 35
1
these reactions are of practical importance, and they could be the rate-determining step for
2
kentonization and aldol condensation in biomass conversion and even decomposition of
3
molecular explosives.73-74 Herein we revisit of the tautomerization of glycine in water.
4
Glycine is the smallest amino acid that has been found in its neutral form (N) in the gas phase
5
or vacuum, while in aqueous solution the zwitterionic form (Z) dominates due to the
6
electrostatic interaction with the water molecules. Although the free energy change from the
7
neutral form to the zwitterionic form is about -7.3 kcal/mol measured in experiments,61, 75 there
8
are two possible reaction paths: the first path is the so called nonassisted concerted mechanism
9
in which the hydroxyl hydrogen (H6) forms a 5-membered ring transition state 1 (TS1) before
10
moving to the amino nitrogen (N1) of the zwitterionic form (Z) (see Fig. 3a). The second path
11
is the assisted concerted mechanism in which the glycine forms a 7-membered ring transition
12
state 2 (TS2) with water and then dehydrates to form a zwitterion structure (Z) (see Fig. 3b).
a)
H6 N1
O4 C3
C2
O5
N
TS1
Z
TS2
Z
Ow
b) Hw H6 N1
O4 C3
C2
O5
13
N
14
Figure 3. Schematic representations and atomic numbering for the neutral, transition, and
15
zwitterionic structure of glycine with: a) nonassisted concerted mechanism and b) assisted
16
concerted mechanism. 12
ACS Paragon Plus Environment
Page 13 of 35 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 2
The free energy profile for the Glycine tautomerization reaction in aqueous solution is
3
depicted in Fig. 4. In the initial state, the grand potential of the system is composed the intrinsic
4
free energy of the quantum system and grand potential energy of the surrounding solvent
5
system, with the external potential corresponding to the molecular geometry and charge
6
distribution in the quantum system. Similar composition holds regarding the grand potentials
7
for the transition state system and for the final state. In each state, the molecular structure and
8
charge distribution of glycine are determined upon quantum DFT calculation, which hereby
9
provide the necessary input information, through the solute-solvent interaction potential, for
10
the CDFT calculation of the hydration free energies.
11 12
Figure 4. Schematic free energy profile for glycine tautomerization reaction in aqueous
13
solution.
14 15
III.
16
1. Intrinsic free energy calculation with QDFT
Calculation Methods and Details
13
ACS Paragon Plus Environment
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 14 of 35
1
To calculate intrinsic reaction free energy ( E R ), we first optimize the structure of
2
glycine in water by treating the solution with polarizable continuum model, and then the
3
obtained geometric configuration is utilized for a single-point frequency calculation in gas
4
phase for determining E R . The quantum calculations are carried out by using the
5
GAUSSIAN 09 package of programs76 at the MP2,77 B3LYP78 and M06-2X79 levels,
6
respectively. Following previous studies of glycine,57, 64, 69 the 6-31+G** basis set has been
7
chosen, and this basis set includes the diffuse functions on all atoms except the hydrogens
8
and the polarization functions on all atoms. Molecular geometries of the reactants, products
9
and transition states (TS1 and TS2) have been optimized using Berny’s algorithms.80 A
10
single-point frequency calculation is followed to ensure that the final structure is obtained
11
without imaginary frequency (reactant and product) or with only a single imaginary
12
frequency (transition states). These TS1 and TS2 transition states are verified by performing a
13
normal-mode analysis, which provides a single imaginary frequency ( υTS1 = -1179 cm-1; υTS2
14
= -1313 cm-1 for the glycine at MP2 level). Intrinsic reaction coordinate (IRC)81 calculations
15
are performed for ensuring the transition states connect to the correct reactants and products.
16
2. Hydration free energy calculation with CDFT
17
2.1 excess free-energy functional
18 19
Following our previous work39, the formulation of the excess free energy functional for inhomogeneous water is formulated as:
20
1 ∆F exc [ ρ (x)]= − k BT ∫∫ ∆ρ ( x1 ) c(2) (x1 , x 2 ;[ ρ0 ])∆ρ ( x 2 ) dx1dx 2 + F B [ ρ ( x )] , 2
21
where ∆ρ (x) = ρ (x) − ρ 0 . On the right hand side of eq.(6), the first term represents the
22
homogeneous reference fluid (HRF) approximation. It amounts to a second-order Taylor
23
expansion of the excess free-energy functional for inhomogeneous water around the bulk one,
14
ACS Paragon Plus Environment
(6)
Page 15 of 35 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
which is known to be equivalent to the hypernetted-chain (HNC) approximation in integral
2
(2) equation theories.22 c (x1 , x2 ;[ρ0 ]) is the two-body direct correlation function (DCF) of the
3
bulk water under the same thermodynamic condition, and it can be extracted from molecular
4
dynamics simulation.82 In this work, the DCF of bulk water serves as an input of the density
5
functional, prepared elaborately in prior.83 The second term in eq.(6) represents the
6
contribution from the multi-body correlations, conventionally designated as bridge functional,
7
which can be rigorously expressed as a series of expansions in terms of the three-body and
8
higher-order correlations of the bulk water. As in many work,39, 82, 84 we here approximate the
9
bridge functional with the one for a reference inhomogeneous hard-sphere fluid system with
10
density distribution ρ (r) :
F B [ ρ ( x )] ≈ FhsB [ ρ ( r )]
11
(7)
12
Here ρ (r ) = ∫ ρ ( x)d Ω . Similar to earlier work,85 the bridge functional of the corresponding
13
hard sphere system is calculated from the modified fundamental measure theory (MFMT)86: exc exc FhsB [ ρ ( r )]=Fhsexc [ ρ ( r )] − Fhs-b ( ρ b ) − µ hs-b ∫ ∆ρ ( r ) d r
14
1 (2) + k BT ∫∫ ∆ρ ( r1 ) chs-b (| r1 − r2 |) ∆ρ ( r2 ) dr1dr2 2
(8)
15
where Fhsexc [ ρ ( r )] is the excess free-energy functional of the inhomogeneous hard-sphere
16
exc exc (2) system; Fhs-b , µ hs-b and chs-b are the excess free-energy, excess chemical potential and
17
(2) DCF for the corresponding hard-sphere bulk system. In this work, Fhsexc [ ρ ( r )] and chs-b are
18
exc exc calculated from MFMT,86 while Fhs-b and µ hs-b are calculated from the Carnahan-Starling
19
equation of state.87 The effective hard-sphere diameter of the reference system is 2.95Å for
20
SPC/E water.33
21
2.2 Numerical aspects
22
The minimization of eq.(5) gives rise to the Euler-Lagrange equation: 15
ACS Paragon Plus Environment
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
1
ρ (x) δ∆F exc [ ρ ( x)] k BT ln = 0. + Vext ( x) + δρ ( x) ρ0
Page 16 of 35
(9)
2
The inhomogeneous density ρ (x) is projected onto an orthorhombic position grid of
3
N x × N y × N z nodes with periodic boundary conditions applied in all three directions. The
4
variational density ρ (x) , subject to the minimization of F [ ρ (x);Vext (x)] , is optimized by the
5
limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method.88 This quasi-Newton
6
algorithm requires the formulation of the functional F [ ρ (x);Vext (x)] and its first derivative
7
with respect to the density at each node. The convolutions in F exc [ ρ (x)] are calculated in
8
k-space by using fast Fourier transforms (FFT) as implemented with the FFTW3 library.89
9
The temperature T = 298.15 K, and water bulk density 0.0333 particles/Å3, corresponding to
10
the mass concentration 0.997 g/cm3.
11
The external potential Vext (x) , involving atomic Lennard-Jones (LJ) and partial charges
12
parameters, is computed only once before the numerical minimization. In each calculation, a
13
single glycine molecule (i.e. in neutral, zwitterionic or transition state form) is fixed at the
14
center of a cubic grid box with box length of approximately 40 Å, and the external potential
15
to the water molecule locating at r with orientation Ω sums up the pair interactions
16
between the site i in glycine and site j in water, reading:
17
σ Vext ( r, Ω ) = ∑∑ 4ε ij ij rij i j a
b
12 6 σ ij qi q j . − + rij 4πε 0 rij
(10)
18
Here rij is the interatomic distance between site i and site j , while qi and q j are the
19
partial charges on these sites. The conventional Lorentz-Berthelot combination rule is used to
20
generate the LJ parameters ε ij and σ ij for the unlike pairs.
16
ACS Paragon Plus Environment
Page 17 of 35 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
For water molecule, we have used the SPC/E model potential parameters.90 The charge
2
distributions in glycine are calculated by fitting the individual (atomic) charges to the natural
3
population analysis (NPA) or molecular electrostatic potential. The latter are computed from
4
the QM wave function at series of points selected by the charges from electrostatic potential
5
(CHELPG) algorithm.91 The LJ parameters of the nonassisted and assisted mechanisms for
6
neutral, zwitterionic and transition state glycine are taken from the OPLS force field database
7
(Tables S1-S6).92
8
IV.
9
1. Molecular Geometric Configurations
Results and Discussion
10
For the nonassisted concerted mechanism, the neutral (N), zwitterionic (Z) and transition
11
state 1 (TS1) conformations of glycine in water are displayed in Fig.3a. In general, the neutral
12
glycine presents a large number of possible conformers in various environments,62, 93-95 and its
13
conformer in water is determined by minimizing the absolute value of its intrinsic free energy
14
by treating water as polarizable continuum.15 Note that the obtained structure of glycine in
15
water is different to its counterpart in gas phase. The same method is applied to determine the
16
structures of glycine in zwitterionic and TS1 states in water. In gas phase, the geometric
17
configurations of zwitterionic and TS1 glycine molecules can’t be located through the
18
minimization of intrinsic free energy,63 and thus the neutral, zwitterionic and transition state 1
19
conformations of glycine are assumed to be identical with their counterparts in water.
20
Fig. 3b presents the configuration involving only one water molecule in the assisted
21
concerted mechanism. The hydroxyl group is oriented to the amino group, which facilitates the
22
formation of intramolecular hydrogen bonds by promoting the proton transfer of H6 atom to
23
Ow atom (assisted mechanism) or N1 atom (nonassisted mechanism). The geometric
24
configurations of glycine in neutral, zwitterionic and transition state 2 (TS2) in the assisted 17
ACS Paragon Plus Environment
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
1
concerted mechanism are determined with the similar method as in the nonassisted concerted
2
mechanism.
3
The geometric parameters described for the neutral (N), zwitterionic (Z) and transition
4
state (TS1 and TS2) of glycine are given in Table 1. All these optimized structures are
5
calculated in the same computational level mentioned in section III. These results are
6
consistent with the observations by other workers for amino acids in water.64-67
7
18
ACS Paragon Plus Environment
Page 18 of 35
Page 19 of 35 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
The Journal of Physical Chemistry
Table 1. Geometric parameters for glycine at the different computational levels used.
N C2-N1 C2-C3 C3-O5 C3-O4 C2-H7 O4-H6 N1-H6 O4-N1 H6-Ow Ow-Hw Hw-N1 N1-C2-C3 O5-C3-C2 O4-C3-C2 N1-H6-O4 O4-H6-Ow H6-Ow-Hw Ow-Hw-N1 N1-C2-C3-O4 C3-C2-N1-H6 C3-C2-N1-Hw
TS1
MP2
B3LYP
M06-2X
MP2
1.46 1.53 1.23 1.34 1.09 1.00 2.57 109.8 123.6 113.8 128.8 0.0 0.0 -
1.47 1.53 1.22 1.33 1.09 1.00 2.57 109.8 123.9 113.4 128.8 0.0 0.0 -
1.46 1.53 1.21 1.33 1.09 0.99 2.57 109.9 123.5 113.9 127.3 0.0 0.0 -
1.48 1.54 1.24 1.31 1.09 1.24 1.27 2.37 103.8 121.8 110.6 140.8 0.0 0.0 -
TS2
B3LYP M06-2X
1.49 1.55 1.23 1.30 1.09 1.26 1.27 2.37 103.8 121.8 110.6 140.0 0.0 0.0 -
1.48 1.55 1.23 1.29 1.09 1.30 1.22 2.37 103.5 120.9 110.9 139.3 0.0 0.0 -
MP2
1.48 1.53 1.24 1.31 1.09 1.17 2.89 1.25 1.24 1.25 114.1 117.6 118.9 170.3 86.0 162.3 -36.0 52.2
19
ACS Paragon Plus Environment
Z
B3LYP M06-2X
1.49 1.54 1.23 1.30 1.09 1.19 2.90 1.22 1.24 1.26 115.2 117.1 119.4 169.8 86.7 162.3 -32.2 47.8
1.48 1.54 1.23 1.29 1.09 1.24 2.87 1.17 1.22 1.27 114.4 117 118.8 168.1 87.9 161.2 -32.6 51.3
MP2
1.50 1.55 1.26 1.27 1.09 1.04 2.56 107.7 115.5 115.1 123.3 -0.5 0.7 -
B3LYP M06-2X
1.50 1.55 1.25 1.27 1.09 1.05 2.55 107.7 116.0 114.9 124.2 -0.3 0.4 -
1.49 1.55 1.24 1.26 1.09 1.05 2.55 107.7 115.6 115.0 122.7 0.0 0.0 -
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
1
In relation to the TS1 structure on the way from neutral to zwitterionic, O4-H6 distance
2
is much larger at M06-2X level (0.31 Å) than at the MP2 (0.24 Å) or B3LYP (0.26 Å). In all
3
cases, the process puts the H6 at a distance of about 1.2 Å of the N1 by extending the O4-H6
4
bond. Moreover, we find that the TS1 glycine is a nearly planar 5-membered ring with very
5
similar distances and bond angles.
6
When a water molecule is involved as part of the reaction mechanism, we find that the
7
Hw links to N1, and the H6 links to Ow. We also observe in Table 1 that the length of the
8
Ow-Hw bond increases to 1.24 Å (MP2 and B3LYP) or 1.22 Å (M06-2X), with the Ow-H6
9
bond length at 1.25 Å (MP2), 1.22 Å (B3LYP) and 1.17 Å (M06-2X), and the O4-H6
10
distance is slightly increased, which allows to release the H6-Ow-Hw water molecule. In the
11
TS2 structure, the formation of a 7-membered nonplanar ring renders the structure less
12
stressed than the TS1 structure, presenting a larger bond angle especially the O4-H6-Ow and
13
Ow-Hw-N1 angles involved in the proton transfer.
14
2. Energetics and Reaction Path
15
Table 2 and Fig. 5 present the relative energies and free energy profiles for glycine at
16
the neutral, zwitterionic, and transition states in water, with different levels of QDFT
17
calculations. In the gas phase for both the assisted and nonassisted mechanisms, the
18
tautomerization process N→Z is clearly unfavorable, with higher activation barriers and very
19
endothermic reaction energies. In aqueous solution, there is a significant reduction of the
20
activation barrier, generating depressed reaction energies of the process being obtained (they
21
even are exothermic for both mechanisms). This is because the dipole moment of glycine
22
increases along the reaction coordinate,72 leading to a continuous decrease of the hydration free
23
energy, and this suggests that the solvent plays an important role in stabilizing the final state.
24
Furthermore, the reduction of the activation barriers in the water-assisted mechanism is
25
slightly more pronounced than that in the non-assisted mechanism, indicating that the
20 ACS Paragon Plus Environment
Page 20 of 35
Page 21 of 35 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
solvent’s role is more significant in the water-assisted mechanism. Indeed, the dipole moment
2
of TS2 is 8.4 D, being larger than 8.2 D for TS1 at the B3LYP level. Tolosa et al.57 proposed a
3
scheme to describe the chemical reaction in water by means of molecular dynamics (MD)
4
simulation combining with quantum calculation at MP2 level, and the MD/MP2 method was
5
applied to study the tautomerization of four amino acids in water using free energy curves at
6
298 K. In their calculation, the TIP3P water model was employed.96
7
The values of the activation and reaction free energies with the nonassisted mechanism
8
and assisted mechanism predicted by RxDFT are quantitatively comparable with the
9
calculation results by Tolosa et al.57 Moreover, the combination of B3LYP with MDFT in the
10
nonassisted mechanism gives very satisfactory predictions (7.83 kcal/mol for the activation
11
free energy and -4.84 kcal/mol for the reaction free energy), compared with Tolosa’s results
12
(6.64 kcal/mol and -4.91 kcal/mol). Besides, the theoretical predictions are overall consistent
13
with the experimental measurements.68, 75, 61, 97
14
Benchmarked with the results from Tolosa, the free energy profiles of glycine in gas
15
phase can be satisfactorily predicted with these three methods, rationalizing our treatments on
16
the geometric structures of glycine in gas phase. In the aqueous solution, the results from
17
B3LYP + MDFT and M06-2X + MDFT method are better than that from MP2 + MDFT
18
method for the nonassisted mechanism, while the free energy profile from MP2 + MDFT
19
method is better than those from B3LYP + MDFT and M06-2X + MDFT methods for the
20
assisted mechanism.
21 22
Table 2. Relative energies (in kcal/mol) for the structures of glycine studied at the different
23
computational levels used. N nonassisted
gas
MP257
TS
0.00 11.56
21 ACS Paragon Plus Environment
Z 20.65
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 22 of 35
MP2
0.00
9.04
20.47
B3LYP
0.00
8.61
20.73
M06-2X
0.00 10.17
23.53
MD/MP257
0.00
6.64
-4.91
MP2 + MDFT
0.00
6.24
-11.38
B3LYP + MDFT
0.00
7.83
-4.84
M06-2X + MDFT
0.00
9.62
-3.08
MP257
0.00 13.40
13.23
MP2
0.00 11.26
11.93
B3LYP
0.00
9.25
14.00
M06-2X
0.00
9.65
15.39
MD/MP257
0.00
5.29
-4.27
MP2 + MDFT
0.00
4.04
-5.85
B3LYP + MDFT
0.00
3.93
-9.08
M06-2X + MDFT
0.00
3.18
-7.22
0.00
6.90
-7.30
solution
gas
assisted
solution
Experimental 1
2 3
Figure 5. Free energy profiles for glycine using a) a nonassisted mechanism, b) an assisted
22 ACS Paragon Plus Environment
Page 23 of 35 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
mechanism. The label G denotes the reaction in gas phase and S denotes the reaction in
2
aqueous solution. The free energy profiles of MD/MP2 method are taken from Ref. 57.
3 4
V.
Conclusions
5
We have proposed an RxDFT by combining QDFT and CDFT, which enable us to
6
address the solvent effect of chemical reactions explicitly. For demonstration, the RxDFT is
7
applied to investigate the reaction paths associated with the tautomerization reaction of
8
glycine in water. In particular, the intrinsic reaction free energy and solvent effect are
9
decoupled, and thereafter treated separately by using QDFT and CDFT. The geometric
10
configurations and intrinsic free energies of and the charge distributions in glycine at
11
different reaction states have been computed with QDFT at the MP2/6-31+G**,
12
B3LYP/6-31+G** and M06-2x/6-31+G** computational level. The obtained geometries and
13
charge distributions determine the glycine-water interaction, which serves, hereby, as an input
14
for the CDFT calculation for determining the change in solvation free energy during the
15
occurrence of tautomerization. The summation of the intrinsic reaction free energy with the
16
solvation free energy difference gives the total reaction free energy in water.
17
The intrinsic free energies are obtained at the MP2, B3LYP and M06-2x computational
18
levels, and the profile confirms that the chemical reaction in gas phase is endothermic. The
19
final free-energy curves for both nonassisted concerted mechanism and assisted concerted
20
mechanism are obtained by using RxDFT, showing that the reaction in water is exothermic.
21
The theoretical predictions from RxDFT are compared with the reported MD/MP2
22
calculations and experimental measurement, and an overall satisfactory agreement has been
23
found, demonstrating that the proposed RxDFT is accurate and reliable.
24
Although water is uniquely considered as a common representative for reaction
25
solvents in this work, the proposed protocol in RxDFT can be extended to address the
23 ACS Paragon Plus Environment
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
1
chemical reactions in other solvents. Due to the low computational cost, we consider RxDFT
2
may provide an efficient alternative to ab initio correlated methods and QM/MM methods for
3
accurate theoretical study of reaction free energies in solvents.
4 5
Associated Content
6
Supporting information
7
The atomic (i.e. site) coordinates, Lennard-Jones parameters and charges for the glycine
8
molecules used in the RxDFT calculation.
9 10
Acknowledgements
11
This work is supported by National Natural Science Foundation of China (No. 91434110 and
12
21878078), National Natural Science Foundation of China for Innovative Research Groups
13
(No. 51621002), and the 111 Project of China (No. B08021). SZ acknowledges the support of
14
Fok Ying Tong Education Foundation (151069).
15 16
References
17
1.
18 19
24. 2.
20 21
3.
Tomasi, J.; Mennucci, B.; Cammi, R., Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999-3093.
4.
24 25
Mennucci, B.; Cammi, R., Continuum Solvation Models in Chemical Physics: From Theory to Applications; John Wiley & Sons, 2007.
22 23
Hirata, F., Molecular Theory of Solvation; Springer Science & Business Media, 2003; Vol.
Reichardt, C.; Welton, T., Solvents and Solvent Effects in Organic Chemistry; John Wiley & Sons, 2011.
5.
Carey, F. A., Organic Chemistry; The McGraw-Hill Companies, 2004.
24 ACS Paragon Plus Environment
Page 24 of 35
Page 25 of 35 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
6.
2 3
Miertus,̃ S.; Tomasi, J., Approximate Evaluations of the Electrostatic Free Energy and Internal Energy Changes in Solution Processes. Chem. Phys. 1982, 65, 239-245.
7.
Miertuš, S.; Scrocco, E.; Tomasi, J., Electrostatic Interaction of a Solute with a Continuum.
4
A Direct Utilizaion of Ab Initio Molecular Potentials for the Prevision of Solvent Effects.
5
Chem. Phys. 1981, 55, 117-129.
6
8.
Prototype SN2 Reaction in Solution: Cl− Attack on CH3Cl. J. Chem. Phys. 2014, 140, 289.
7 8 9
Kuechler, E. R.; York, D. M., Quantum Mechanical Study of Solvent Effects in a
9.
Car,
R.;
Parrinello,
M.,
Unified
Approach
for
Molecular
Dynamics
and
Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471.
10
10. Vreven, T.; Frisch, M. J.; Kudin, K. N.; Schlegel, H. B.; Morokuma, K., Geometry
11
Optimization with QM/MM Methods II: Explicit Quadratic Coupling. Mol. Phys. 2006,
12
104, 701-714.
13 14 15 16
11. Teodoro, L.; Fawzi, M.; Alessandro, L.; Parrinello, M., An Efficient Real Space Multigrid QM/MM Electrostatic Coupling. J. Chem. Theory Comput. 2005, 1, 1176-1184. 12. Orozco, M.; Luque, F. J., Theoretical Methods for the Description of the Solvent Effect in Biomolecular Systems. Chem. Rev. 2000, 100, 4187-4226.
17
13. Muñoz, J.; Barril, X.; Luque, F. J.; Gelpí, J. L.; Orozco, M., Partitioning of Free Energies
18
of Solvation into Fragment Contributions: Applications in Drug Design; Springer US,
19
2001, p 143-168.
20 21
14. Vreven, T.; Morokuma, K., On the Application of the IMOMO (Integrated Molecular Orbital + Molecular Orbital) Method. J. Comput. Chem. 2015, 21, 1419-1432.
22
15. Vreven, T.; Byun, K. S.; Komáromi, I.; Dapprich, S.; Montgomery, J. A.; Morokuma, K.;
23
Frisch, M. J., Combining Quantum Mechanics Methods with Molecular Mechanics
24
Methods in ONIOM. J. Chem. Theory Comput. 2006, 2, 815-826.
25 ACS Paragon Plus Environment
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
1
16. Mikkelsen, K. V.; Jo/Rgensen, P.; Jensen, H. J. R. A., A Multiconfiguration
2
Self-Consistent Reaction Field Response Method. J. Chem. Phys. 1994, 100, 6597-6607.
3
17. Bishop, D. M., Aspects of Non-Linear-Optical Calculations. Adv. Quantum Chem. 1994,
4 5 6
25, 1-45. 18. Shelton, D. P.; Rice, J. E., Measurements and Calculations of the Hyperpolarizabilities of Atoms and Small Molecules in the Gas Phase. Chem. Rev. 1994, 94, 3-29.
7
19. Galván, I. F.; Sánchez, M. L.; Martı́N, M. E.; Valle, F. J. O. D.; Aguilar, M. A., ASEP/MD:
8
A Program for the Calculation of Solvent Effects Combining QM/MM Methods and the
9
Mean Field Approximation. Comput. Phys. Commun. 2003, 155, 244-259.
10
20. Sanchez, M. L.; Aguilar, M. A.; Olivares, d. V., F. J, Study of Solvent Effects by Means of
11
Averaged Solvent Electrostatic Potentials Obtained from Molecular Dynamics Data. J.
12
Comput. Chem. 2015, 18, 313-322.
13
21. Sánchez, M. L.; Martín, M. E.; Aguilar, M. A.; Olivares, d. V., F. J, Solvent Effects by
14
Means of Averaged Solvent Electrostatic Potentials: Coupled Method. J. Comput. Chem.
15
2015, 21, 705-715.
16
22. Hansen, J. P.; Mcdonald, I. R., Theory of Simple Liquids (Second Edition), 1990.
17
23. Hirata, F.; Sato, H.; Ten-No, S.; Kato, S., The RISM-SCF/MCSCF Approach for
18
Chemical Processes in Solutions. Eastern Hemisphere Distribution 2001, 417.
19
24. Andersen, H. C.; Chandler, D., Optimized Cluster Expansions for Classical Fluids. I.
20
General Theory and Variational Formulation of the Mean Spherical Model and Hard
21
Sphere Percus-Yevick Equations. J. Chem. Phys. 1972, 57, 1918-1929.
22 23
25. Chandler, D.; Andersen, H. C., Optimized Cluster Expansions for Classical Fluids. II. Theory of Molecular Liquids. J. Chem. Phys. 1972, 57, 1930-1937.
26 ACS Paragon Plus Environment
Page 26 of 35
Page 27 of 35 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
26. Andersen, H. C.; Chandler, D.; Weeks, J. D., Optimized Cluster Expansions for Classical
2
Fluids. III. Applications to Ionic Solutions and Simple Liquids. J. Chem. Phys. 1972, 57,
3
2626-2631.
4 5
27. Hudson, S.; Andersen, H. C., Optimized Cluster Expansions for Classical Fluids. IV. Primitive Model Electrolyte Solutions. J. Chem. Phys. 1974, 60, 2188-2188.
6
28. Wesolowski, T. A.; Weber, J., Kohn-Sham Equations with Constrained Electron Density:
7
An Iterative Evaluation of the Ground-State Electron Density of Interacting Molecules.
8
Chem. Phys. Lett. 1996, 248, 71-76.
9
29. Naka, K.; Sato, H.; Morita, A.; Hirata, F.; Kato, S., Rism-Scf Study of the Free-Energy
10
Profile of the Menshutkin-Type Reaction NH3 + CH3Cl → NH3CH3+ + Cl− in Aqueous
11
Solution. Theor. Chem. Acc. 1999, 102, 165-169.
12
30. Zhao, S.; Liu, Y.; Chen, X.; Lu, Y.; Liu, H.; Hu, Y., Chapter One-Unified Framework of
13
Multiscale Density Functional Theories and Its Recent Applications. Advances in
14
Chemical Engineering 2015, 47, 1-83.
15 16 17 18
31. Wu, J. Z., Density Functional Theory for Liquid Structure and Thermodynamics; Springer Berlin Heidelberg, 2009. 32. Zhao, S. L.; Liu, Y.; Liu, H. L.; Wu, J. Z., Site-Site Direct Correlation Functions for Three Popular Molecular Models of Liquid Water. J. Chem. Phys. 2013, 139, 064509.
19
33. Liu, Y.; Zhao, S.; Wu, J., A Site Density Functional Theory for Water: Application to
20
Solvation of Amino Acid Side Chains. J. Chem. Theory Comput. 2013, 9, 1896-1908.
21
34. Lian, C.; Cai, C.; Shen, X.; Zhao, S.; Yu, X.; Liu, H., Improved Oxidation of Hydrogen
22
Off-Gas by Hydrophobic Surface Modification: A Multiscale Density Functional Theory
23
Study. Particuology 2018.
24 25
35. Wu, J., Classical Density Functional Theory for Molecular Systems; Springer Singapore, 2017.
27 ACS Paragon Plus Environment
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
1
36. Chandler, D.; McCoy, J. D.; Singer, S. J., Density Functional Theory of Nonuniform
2
Polyatomic Systems. I. General Formulation. J. Chem. Phys. 1986, 85, 5971-5976.
3
37. Reddy, G.; Lawrence, C. P.; Skinner, J. L.; Yethiraj, A., Liquid State Theories for the
4
Structure of Water. J. Chem. Phys. 2003, 119, 13012-13016.
5
38. Sundararaman, R.; Arias, T. A., Efficient Classical Density-Functional Theories of
6
Rigid-Molecular Fluids and a Simplified Free Energy Functional for Liquid Water.
7
Comput. Phys. Commun. 2014, 185, 818-825.
8 9 10 11
39. Zhao, S.; Jin, Z.; Wu, J., New Theoretical Method for Rapid Prediction of Solvation Free Energy in Water. J. Phys. Chem. B 2011, 115, 6971-6975. 40. Wu, H.; Li, Y.; Kadirov, D.; Zhao, S.; Lu, X.; Liu, H., An Efficient Molecular Approach for Quantifying Solvent-Mediated Interaction. Langmuir 2017, 33, 11817-11824.
12
41. Liu, Y.; Liu, H.; Hu, Y.; Jiang, J., Development of a Density Functional Theory in
13
Three-Dimensional Nanoconfined Space: H2 Storage in Metal-Organic Frameworks. J.
14
Phys. Chem. B 2009, 113, 12326-12331.
15 16
42. Liu, Y.; Liu, H.; Hu, Y.; Jiang, J., Density Functional Theory for Adsorption of Gas Mixtures in Metal-Organic Frameworks. J. Phys. Chem. B 2010, 114, 2820.
17
43. Cheng, L.; Zhao, S.; Liu, H.; Wu, J., Time-Dependent Density Functional Theory for the
18
Charging Kinetics of Electric Double Layer Containing Room-Temperature Ionic Liquids.
19
J. Chem. Phys. 2016, 145, 204707.
20 21 22 23
44. Zhao, S.; Wu, J., Self-Consistent Equations Governing the Dynamics of Nonequilibrium Colloidal Systems. J. Chem. Phys. 2011, 134, 054514. 45. Ma, M.; Zhao, S.; Liu, H.; Xu, Z., Microscopic Insights into the Efficiency of Capacitive Mixing Process. AlChE J. 2017, 63, 1785-1791.
28 ACS Paragon Plus Environment
Page 28 of 35
Page 29 of 35 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
46. Ma, M.; Zhao, S.; Xu, Z., Investigation of Dielectric Decrement and Correlation Effects
2
on Electric Double-Layer Capacitance by Self-Consistent Field Model. Communications
3
in Computational Physics 2016, 20, 441-458.
4 5 6 7 8 9
47. Petrosyan, S. A.; Briere, J. F.; Roundy, D.; Arias, T. A., Joint Density-Functional Theory for Electronic Structure of Solvated Systems. Physical Review B 2007, 75. 48. Sundararaman, R. Joint Density-Functional Methods for First-Principles Chemistry in Solution. Cornell University, 2013. 49. Letchworth-Weaver, K.; Sundararaman, R.; A. Arias, T., First Principles Free-Energy Theory of Solvation with Atomic Scale Liquid Structure, 2017.
10
50. Sundararaman, R.; Letchworth-Weaver, K.; Schwarz, K. A.; Gunceler, D.; Ozhabes, Y.;
11
Arias, T. A., JDFTx: Software for Joint Density-Functional Theory. SoftwareX 2017, 6,
12
278-284.
13
51. Deng, S.; Cai, W., Extending the Fast Multipole Method for Charges inside a Dielectric
14
Sphere in an Ionic Solvent: High Order Image Approximations for Reaction Fields.
15
Journal of Computational Physics 2007, 227, 1246-1266.
16
52. Lin, Y.; Baumketner, A.; Deng, S.; Xu, Z.; Jacobs, D.; Cai, W., An Image-Based Reaction
17
Field Method for Electrostatic Interactions in Molecular Dynamics Simulations of
18
Aqueous Solutions. J. Chem. Phys. 2009, 131, 10B608.
19 20
53. Lin, Y.; Baumketner, A.; Song, W.; Deng, S.; Jacobs, D.; Cai, W., Ionic Solvation Studied by Image-Charge Reaction Field Method. J. Chem. Phys. 2011, 134, 01B627.
21
54. Qin, P.; Xu, Z.; Cai, W.; Jacobs, D., Image Charge Methods for a Three-Dielectric-Layer
22
Hybrid Solvation Model of Biomolecules. Communications in Computational Physics
23
2009, 6, 955-977.
29 ACS Paragon Plus Environment
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
1
55. Berg, M.; Harris, A. L.; Harris, C. B., Rapid Solvent-Induced Recombination and Slow
2
Energy Relaxation in a Simple Chemical Reaction: Picosecond Studies of Iodine
3
Photodissociation in CCl4. Phys. Rev. Lett. 1985, 54, 951-954.
4
56. Rosenthal, S. J.; Xie, X.; Du, M.; Fleming, G. R., Femtosecond Solvation Dynamics in
5
Acetonitrile: Observation of the Inertial Contribution to the Solvent Response. J. Chem.
6
Phys. 1991, 95, 4715-4718.
7
57. Tolosa, S.; Hidalgo, A.; Sansón, J. A., Amino Acid Tautomerization Reactions in
8
Aqueous Solution Via Concerted and Assisted Mechanisms Using Free Energy Curves
9
from Md Simulation. J. Phys. Chem. B 2012, 116, 13033-13044.
10 11 12 13 14 15
58. Zhao, S.; Ramirez, R.; Vuilleumier, R.; Borgis, D., Molecular Density Functional Theory of Solvation: From Polar Solvents to Water. J. Chem. Phys. 2011, 134, 194102. 59. Laidler, K. J.; King, M. C., Development of Transition-State Theory. J. Phys. Chem. 1983, 87, 2657-2664. 60. Smith, M. B.; March, J., March's Advanced Organic Chemistry: Reactions, Mechanisms, and Structure; John Wiley & Sons, 2007.
16
61. Wada, G.; Tamura, E.; Okina, M.; Nakamura, M., On the Ratio of Zwitterion Form to
17
Uncharged Form of Glycine at Equilibrium in Various Aqueous Media. Bull. Chem. Soc.
18
Jpn. 1982, 55, 3064-3067.
19 20
62. Császár, A. G., On the Structures of Free Glycine and α-Alanine. J. Mol. Struct. 1995, 346, 141-152.
21
63. Tortonda, F. R.; Pascual-Ahuir, J. L.; Silla, E.; Tuñón, I., Why Is Glycine a Zwitterion in
22
Aqueous Solution? A Theoretical Study of Solvent Stabilising Factors. Chem. Phys. Lett.
23
1996, 260, 21-26.
30 ACS Paragon Plus Environment
Page 30 of 35
Page 31 of 35 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
64. Tortonda, F. R.; Pascualahuir, J. L.; Silla, E.; Tuñón, I.; Ramı́Rez, F. J., Aminoacid
2
Zwitterions in Solution: Geometric, Energetic, and Vibrational Analysis Using Density
3
Functional Theory-Continuum Model Calculations. J. Chem. Phys. 1998, 109, 592-603.
4
65. Balta, B.; Aviyente, V., Solvent Effects on Glycine. I. A Supermolecule Modeling of
5
Tautomerization Via Intramolecular Proton Transfer. J. Comput. Chem. 2003, 24,
6
1789-802.
7
66. Tortonda, F. R.; Silla, E.; Tuñón, I.; Rinaldi, D.; Ruiz-López, M. F., Intramolecular Proton
8
Transfer of Serine in Aqueous Solution. Mechanism and Energetics. Theor. Chem. Acc.
9
2000, 104, 89-95.
10
67. Bowman, J. C.; Shadwick, B. A.; Morrison, P. J., Theoretical Study of the
11
Tautomeric/Conformational Equilibrium of Aspartic Acid Zwitterions in Aqueous
12
Solution. J. Phys. Chem. A 2000, 104, 6834-6843.
13 14 15 16
68. Slifkin, M.; Ali, S., Thermodynamic Parameters of the Activation of Glycine Zwitterion Protonation Reactions. J. Mol. Liq. 1984, 28, 215-221. 69. Tuñón, I.; Silla, E.; Ruiz-López, M. F., On the Tautomerization Process of Glycine in Aqueous Solution. Chem. Phys. Lett. 2000, 321, 433-437.
17
70. Senn, H. M.; Margl, P. M.; Schmid, R.; Ziegler, T.; Blöchl, P. E., Ab Initio Molecular
18
Dynamics with a Continuum Solvation Model. J. Chem. Phys. 2003, 118, 1089-1100.
19
71. Bandyopadhyay, P.; Gordon, M. S.; Mennucci, B.; Tomasi, J., An Integrated Effective
20
Fragment-Polarizable Continuum Approach to Solvation: Theory and Application to
21
Glycine. J. Chem. Phys. 2002, 116, 5023-5032.
22 23
72. Leung, K.; Rempe, S. B., Ab Initio Molecular Dynamics Study of Glycine Intramolecular Proton Transfer in Water. J. Chem. Phys. 2005, 122, 184506.
31 ACS Paragon Plus Environment
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
1
73. Resasco, D. E.; Wang, B.; Crossley, S., Zeolite-Catalysed C-C Bond Forming Reactions
2
for Biomass Conversion to Fuels and Chemicals. Catalysis Science & Technology 2016, 6,
3
2543-2559.
4
74. Wang, B.; Wright, D.; Cliffel, D.; Haglund, R.; Pantelides, S. T., Ionization-Enhanced
5
Decomposition of 2, 4, 6-Trinitrotoluene (Tnt) Molecules. J. Phys. Chem. A 2011, 115,
6
8142-8146.
7
75. Edsall, J. T.; Blanchard, M. H., The Activity Ratio of Zwitterions and Uncharged
8
Molecules in Ampholyte Solutions. The Dissociation Constants of Amino Acid Esters. J.
9
Am. Chem. Soc 1933, 55, 2337-2353.
10
76. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J.
11
R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A., Gaussian 09 Revision A.1.
12
Gaussian Inc. 2009.
13 14
77. Head-Gordon, M.; Pople, J. A.; Frisch, M. J., MP2 Energy Evaluation by Direct Methods. Chem. Phys. Lett. 1988, 153, 503-506.
15
78. Stephens, P.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J., Ab-Intio Calculation of
16
Vibrational Absorption and Circular-Dichroism Spectra Using Density-Function
17
Force-Fields. 1994.
18
79. Zhao, Y.; Truhlar, D. G., The M06 Suite of Density Functionals for Main Group
19
Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States,
20
and Transition Elements: Two New Functionals and Systematic Testing of Four
21
M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215-241.
22
80. Schlegel, H. B., Optimization of Equilibrium Geometries and Transition Structures. J.
23 24 25
Comput. Chem. 1982, 3, 214-218. 81. Fukui, K., The Path of Chemical Reactions - the IRC Approach. Acc. Chem. Res. 1981, 14, 471-476.
32 ACS Paragon Plus Environment
Page 32 of 35
Page 33 of 35 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 2 3 4
82. Zhao, S.; Wu, J., An Efficient Method for Accurate Evaluation of the Site-Site Direct Correlation Functions of Molecular Fluids. Mol. Phys. 2011, 109, 2553-2564. 83. Zhao, S.; Liu, H.; Ramirez, R.; Borgis, D., Accurate Evaluation of the Angular-Dependent Direct Correlation Function of Water. J. Chem. Phys. 2013, 139, 034503.
5
84. Levesque, M.; Vuilleumier, R.; Borgis, D., Scalar Fundamental Measure Theory for Hard
6
Spheres in Three Dimensions: Application to Hydrophobic Solvation. J. Chem. Phys.
7
2012, 137, 034115.
8 9 10 11 12 13 14 15 16 17
85. Zhao, S.; Jin, Z.; Wu, J., New Theoretical Method for Rapid Prediction of Solvation Free Energy in Water. J. Phys. Chem. B 2011, 115, 6971-6975. 86. Yu,
Y.
X.; Wu,
J.,
Structures
of
Hard-Sphere
Fluids
from a
Modified
Fundamental-Measure Theory. J. Chem. Phys. 2002, 117, 10156-10164. 87. Carnahan, N. F.; Starling, K. E., Equation of State for Nonattracting Rigid Spheres. J. Chem. Phys. 1969, 51, 635-636. 88. Byrd, R. H.; Lu, P.; Nocedal, J.; Zhu, C., A Limited Memory Algorithm for Bound Constrained Optimization. Siam Journal on Scientific Computing 1995, 16, 1190-1208. 89. Frigo, M.; Johnson, S. G., The Design and Implementation of FFTW3. Proceedings of the IEEE 2005, 93, 216-231.
18
90. Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F. V.; Hermans, J., Interaction
19
Models for Water in Relation to Protein Hydration; Springer Netherlands, 1981, p
20
331-342.
21
91. Breneman, C. M.; Wiberg, K. B., Determining Atom-Centered Monopoles from
22
Molecular Electrostatic Potentials. The Need for High Sampling Density in Formamide
23
Conformational Analysis. J. Comput. Chem. 1990, 11, 361-373.
33 ACS Paragon Plus Environment
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
1
92. Jorgensen, W. L.; And, D. S. M.; Tiradorives, J., Development and Testing of the OPLS
2
All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J.
3
Am. Chem. Soc. 1996, 118, 11225-11236.
4 5
93. Barone, V.; Adamo, C.; Lelj, F., Conformational Behavior of Gaseous Glycine by a Density Functional Approach. J. Chem. Phys. 1995, 102, 364-370.
6
94. Nguyen, D. T.; Scheiner, A. C.; Andzelm, J. W.; Sirois, S.; Salahub, D. R.; Hagler, A. T., A
7
Density Functional Study of the Glycine Molecule: Comparison with Post-Hartree-Fock
8
Calculations and Experiment. J. Comput. Chem. 2015, 18, 1609-1631.
9
95. Császár, A. G., Conformers of Gaseous α-Alanine. J. Phys. Chem. 1996, 100, 3541-3551.
10
96. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L.,
11
Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys.
12
1983, 79, 926-935.
13 14
97. Haberfield, P., What Is the Energy Difference between H2NCH2CO2H and +H3NCH2CO2-? J. Chem. Educ. 1980, 57, 346-347.
15
34 ACS Paragon Plus Environment
Page 34 of 35
Page 35 of 35 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 2
TOC Graphic
35 ACS Paragon Plus Environment