ABSTRACT

Title of dissertation: ANALYZING DYNAMICAL PROCESSES
WITH LOCAL MOLECULAR FIELD THEORY

Renjie Zhao, Doctor of Philosophy, 2023

Dissertation directed by: Professor John D. Weeks
Department of Chemistry and Biochemistry
Institute for Physical Science and Technology

Local molecular field (LMF) theory provides a framework for describing the

collective response of a system to long-range interactions in nonuniform liquids.

Based on this theory, different roles played by the short and long-range components

of the intermolecular interactions can be disentangled in determining relevant struc-

tural and thermodynamic properties in equilibrium. Furthermore, in dynamical

processes, nonlocal long-range interactions are often associated with long relaxation

times, and can contribute significantly to the stability of the system in different

phases. In this thesis, LMF theory is utilized to quantify and analyze the dynami-

cal effects arising from long-range Coulomb interactions in aqueous solutions, while

elucidating how they are connected to strong local forces and fluctuations.

The first half of the work concerns ionic and dipolar solvation dynamics, which

plays an essential role in many solution phase chemical reactions. The physical

models of Gaussian-smoothed charge and dipole distributions are conceptualized

from LMF theory to investigate the molecular origins of linear and nonlinear effects



in solvation dynamics. The long-range component of the solute-solvent electrostatic

interaction is shown to underlie the linear response behavior of the system, while

the short-range interactions introduce additional nonlinear effects. The LMF-based

solvation models further demonstrate their functionality in probing the intrinsic

dielectric dispersion of solvent water. The second half of the work is focused on

the nucleation processes in the aqueous environment. Simulating crystal nucleation

from solutions requires efficient treatments for intermolecular interactions to drive

the transitions on time scales affordable to molecular dynamics simulations. For

this purpose, a LMF-based molecular model is employed to capture the renormalized

long-range interactions, and well-tempered metadynamics is adopted to enhance the

fluctuations arising from short-range interactions. By comparing to a short-range

reference model, the necessity of long-range interactions in explaining metastability

is revealed. Temporal fluctuations and direct evidence for the two-step nucleation

mechanism are observed through the analysis using a deep learning-based approach.

The results about these two types of dynamical processes contribute to a deeper

understanding of the roles of short and long-ranges interactions in the aqueous

systems.



ANALYZING DYNAMICAL PROCESSES WITH LOCAL
MOLECULAR FIELD THEORY

by

Renjie Zhao

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy

2023

Advisory Committee:
Professor John D. Weeks, Chair/Advisor
Professor Pratyush Tiwary, Co-chair
Professor Christopher Jarzynski
Professor Garegin Papoian
Professor Yifei Mo, Dean’s Representative



© Copyright by
Renjie Zhao

2023



Acknowledgments

I would like to greatly thank Dr. John Weeks, Dr. Pratyush Tiwary, and Dr.

Rick Remsing for their guidance and help in completing the research projects during

my Ph.D. study.

I would also like to sincerely thank Dr. James W Evans and Dr. Chenglin

Luo for their selfless help in the difficult times of my life.

Many thanks to Dr. Yin li for his encouraging in the good old days and for

sharing his academic taste with me.

I would also like to sincerely thank Xiaofan Wu.

Finally, I would like to express to my parents the gratitude that transcends

any words.

ii



Table of Contents

List of Figures v

List of Abbreviations and Symbols vii

List of Publications from Doctoral Research ix

1 Introduction 1
1.1 Review of Weeks-Chandler-Andersen Theory . . . . . . . . . . . . . . 3
1.2 Review of Local Molecular Field Theory . . . . . . . . . . . . . . . . 6
1.3 Chapter Overview of Thesis . . . . . . . . . . . . . . . . . . . . . . . 10

2 Solvation Dynamics in Water I: Origins of Linear and Nonlinear Dynamic
Responses from the Local Molecular Field Theory Perspective 14
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2 Linear Response to Gaussian Charge Distributions in Equilibrium . . 18
2.3 Linear Response Theory for Solvation Dynamics . . . . . . . . . . . . 22
2.4 Molecular Origins of Nonlinear Effects . . . . . . . . . . . . . . . . . 27
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3 Solvation Dynamics in Water II: Intrinsic Dielectric Dispersion of Water 35
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 Solvation of Dipolar Solutes in Equilibrium . . . . . . . . . . . . . . . 37

3.2.1 Bell Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2.2 Gaussian-Smoothed Dipole Model . . . . . . . . . . . . . . . . 40
3.2.3 Model-Based Dielectric Propensities . . . . . . . . . . . . . . . 47

3.3 Linear Response Description for Relaxation Time Scales in Solvation
Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.3.1 linear Irreversible Processes . . . . . . . . . . . . . . . . . . . 49
3.3.2 Probing the Dielectric Dispersion of Water with Gaussian-

Smoothed Solute Models . . . . . . . . . . . . . . . . . . . . . 55
3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

iii



4 Urea Nucleation in Water I: Exploring Free Energy Landscape with Local
Molecular Field Theory and Well-tempered Metadynamics 60
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 Local Molecular Field Theory-Based Models . . . . . . . . . . . . . . 66

4.2.1 Gaussian-truncated Model . . . . . . . . . . . . . . . . . . . . 66
4.2.2 Short Solvent Model . . . . . . . . . . . . . . . . . . . . . . . 69

4.3 Well-tempered Metadynamics . . . . . . . . . . . . . . . . . . . . . . 72
4.4 Sampling Nucleation Events and Constructing Free Energy Surfaces . 77
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

5 Urea Nucleation in Water II: Analysis for Metastability and Nucleation Mech-
anism 86
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.2 State Predictive Information Bottleneck . . . . . . . . . . . . . . . . 89
5.3 Collective Variables for Characterizing Nucleation of Urea in Water . 93
5.4 Influence of Long-range Interactions on Metastability . . . . . . . . . 94
5.5 Two-step Nucleation Mechanism . . . . . . . . . . . . . . . . . . . . . 97

5.5.1 Temporal Fluctuations . . . . . . . . . . . . . . . . . . . . . . 98
5.5.2 Nucleation Pathways . . . . . . . . . . . . . . . . . . . . . . . 99

5.6 Interpretation of Reaction Coordinates . . . . . . . . . . . . . . . . . 101
5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

6 Conclusions and Future Work 107

A Simulation Details for Solvation Dynamics 110

B Simulation Details for Urea Nucleation in Water 111

Bibliography 113

iv



List of Figures

1.1 Illustration of the force cancellation argument . . . . . . . . . . . . . 4
1.2 Breakdown of the force cancellation argument illustrated in a nonuni-

form system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Separation of the Coulomb potential . . . . . . . . . . . . . . . . . . 9

2.1 Probability distribution of the interaction energy between a Gaussian
charge distribution and bulk water . . . . . . . . . . . . . . . . . . . 22

2.2 Solvation dynamics for exciting a Gaussian charge distribution with
appropriately chosen width . . . . . . . . . . . . . . . . . . . . . . . . 28

2.3 Solvation dynamics for exciting a Gaussian charge distribution with
inadequate width . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.4 Solvation dynamics for exciting a point charge inside a hard solute core 32

3.1 Comparison between the point dipole potential and the the Gaussian
dipole potential in radial direction . . . . . . . . . . . . . . . . . . . . 41

3.2 Charge density of Gaussian dipole distribution . . . . . . . . . . . . . 45
3.3 Charging free energy of a Gaussian dipole with different smoothing

lengths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.4 Linear response description of solvation dynamics for exciting Gaus-

sian charge and Gaussian dipole distributions . . . . . . . . . . . . . 56
3.5 Comparison between linear response theory-based model and Born-

like models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.1 Illustration of free energy change described by CNT . . . . . . . . . . 62
4.2 Oxygen-oxygen pair distribution functions of SPC/E water for the

full model and the GT model . . . . . . . . . . . . . . . . . . . . . . 68
4.3 Effective long-range solute-solute interactions in the SS model for

carbon and oxygen sites of urea in the aqueous solution . . . . . . . . 72
4.4 Illustration of the working principle of metadynamics . . . . . . . . . 75
4.5 Characteristic vectors and corresponding intermolecular angles de-

fined for two neighboring urea molecules . . . . . . . . . . . . . . . . 78
4.6 Time series of averaged intermolecular angles . . . . . . . . . . . . . . 79
4.7 Free energy surfaces in the space of averaged intermolecular angles . . 80

v



4.8 Different structures of urea in the aqueous solution . . . . . . . . . . 81
4.9 Comparison among the free energy differences between the liquid-like

state and form I computed with different models . . . . . . . . . . . . 82

5.1 Illustration of the neural network configuration for SPIB . . . . . . . 92
5.2 Free energy barrier between the liquid-like basin and the primary

crystal basin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
5.3 State labels learnt by SPIB for the SS model . . . . . . . . . . . . . . 95
5.4 Comparison between state labels learnt by SPIB for the GT model

and the full model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.5 Cluster size information during transitions between the liquid-like

basin and the primary crystal basin . . . . . . . . . . . . . . . . . . . 100
5.6 State labels in the SPIB reaction coordinate space for the SS model . 101
5.7 Normalized coefficients for the collective variables in the SPIB-learned

reaction coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.8 Time series of averaged intermolecular angles and pair orientational

entropies in the same trajectory . . . . . . . . . . . . . . . . . . . . . 103

vi



List of Abbreviations and Symbols

CNT Classical Nucleation Theory
CV Collective variable
FES Free energy surface
GT Gaussian-truncated
LAMMPS Large-scale Atomic/Molecular Massively Parallel Simulator
LJ Lennard-Jones
LMF Local Molecular Field
MD Molecular dynamics
PMF Potential of mean force
RC Reaction coordinate
SPC/E Extended simple point charge
SPIB State Predictive Information Bottleneck
TCF Time correlation function
WCA Weeks-Chandler-Andersen
WTmetaD Well-tempered metadynamics
YBG Yvon-Born-Green

β inverse thermal energy 1/kBT
ε dielectric constant or dielectric permittivity
εh high frequency limit of the dielectric permittivity
εLJ Lennard-Jones energy parameter
θ̄ averaged intermolecular angle
µ dipole moment
ρ number density
ρ∗ number density in reduced units
ρq induced solvent charge density
ρQ solute charge density
σ width of the Gaussian distribution
σLJ Lennard-Jones range parameter
τD Debye relaxation time
τL solvent longitudinal relaxation time
χqq charge-charge correlation function
Ω grand potential

C(t) normalized time correlation function
erf error function
erfc complementary error function
f(t) phase space distribution function
g(r) pair correlation function

vii



H Hamiltonian
H0 Hamiltonian of the unperturbed system
H′ perturbation Hamiltonian
kB Boltzmann constant
l width of the Gaussian distribution
N8+ population of molecules with coordination numbers greater than 8
N11+ population of molecules with coordination numbers greater than 11
Rc radius of the hard sphere solute
S(r, θ) pair orientational entropy
S(t) normalized dynamic Stokes shift
T temperature
T ∗ temperature in reduced units
u0(r) short-range repulsive component of the Lennard-Jones potential
u1(r) long-range attractive component of the Lennard-Jones potential
uLJ(r) Lennard-Jones potential
v0(r) short-range component of the Coulomb interaction
v1(r) long-range component of the Coulomb interaction
vQ electrostatic potential due to the solute

viii



List of Publications from Doctoral Research

1. Zhao, R., Remsing, R. C., and Weeks, J. D. (2020). Response theory for static
and dynamic solvation of ionic and dipolar solutes in water. Journal of Statistical
Physics, 180, 721–738.

2. Wang, D., Zhao, R., Weeks, J. D., and Tiwary, P. (2022). Influence of long-range
forces on the transition states and dynamics of NaCl ion-pair dissociation in water.
The Journal of Physical Chemistry B, 126 (2), 545-551.*

3. Zhao, R., Zou, Z., Weeks, J. D., and Tiwary, P. (2023). Quantifying the relevance
of long-range forces for crystal nucleation in water. arXiv preprint arXiv:2307.13189.

*The main content of this paper is not included in this thesis.

ix



Chapter 1: Introduction

Understanding the collective behaviors of complex molecular systems relies

on statistical mechanical treatments with idealized molecular models that capture

the essential structure and energetics of the microscopic constituents of the systems.

Modern experimental techniques have facilitated the acquirement of the microscopic

information of physical systems, but the time and length-scale limitations may still

hinder the direct observation of molecular details in many complex dynamical pro-

cesses. Furthermore, the conceptual deconstruction of molecular models inevitably

requires theoretical and computational approaches. For instance, Local molecular

field (LMF) theory suggests that different roles are played by the short and long-

range components of the intermolecular interactions in nonuniform liquids [1–5]. It

has been demonstrated for relevant equilibrium structural [1–3,5,6] and thermody-

namic [7–9] properties in nonuniform polar liquids through the combined approach of

statistical mechanics and molecular dynamics (MD) simulations. When it comes to

nonequilibrium processes, one can naturally extent the LMF framework to modelling

and analyzing the dynamic behaviors of complex liquid systems. In this thesis, two

kinds of nonequilibrium phenomena, solvation dynamics and crystal nucleation in

the aqueous environment, are investigated with MD simulations and the theoretical

1



tools developed based on LMF theory.

In a liquid system, intermolecular interactions can typically be separated into

short-range interactions (such as the repulsive excluded volume interaction and hy-

drogen bonding) that account for the local structure of the liquid, and long-range

interactions (such as the long-range tail of Coulomb interaction) that are essential

for interpreting dielectric screening and other nonlocal effects. LMF theory provides

a physically suggestive framework for describing the different effects attributed to

short and long-range interactions in nonuniform liquids. In nonequilibrium pro-

cesses, the influences of long-range interactions can manifest in more intricate ways.

For instance, fluctuations arising from long-range interactions that decay at large

distances with a power smaller than the spatial dimension can often exhibit slow

relaxation [10,11]. Therefore, advanced sampling methods are in need for quantify-

ing the effects of long-range interactions that possibly emerge on a long timescale

relative to conventional MD simulations. Moreover, long-range interactions can also

influence the stability of solvated species in the aqueous environment [12], thereby

requiring analysis tools for examining the stability of different structures. These

subtleties will be addressed in this thesis through studying the crystal nucleation of

urea in water.

In the remainder of the introduction, LMF theory and the ideas underlying

this theory are reviewed, followed by an outline of the subsequent chapters of the

thesis.

2



1.1 Review of Weeks-Chandler-Andersen Theory

The physical insight used by LMF theory for treating the average effects of

the long-range interactions in nonuniform fluids can be traced backed to the key

idea which was made by Widom to to explain the van der Waals equation for simple

monatomic liquids [13] and then further quantitatively implemented into the pertur-

bation theory of dense uniform liquids developed by Weeks, Chandler, and Andersen

(WCA) [14, 15]. Lies at the heart of this idea is the force cancellation argument,

which states that the attractive intermolecular forces nearly cancel in dense uniform

liquids.

The scenario described by the force cancellation argument can be illustrated by

considering a dense uniform Lennard-Jones (LJ) fluid away from the critical point.

The interaction between the constituting particles is described by the LJ potential

uLJ(r) = 4εLJ

[(σLJ

r

)12
−
(σLJ

r

)6]
, (1.1)

where εLJ measures the depth of the potential well and σLJ is the zero-crossing

distance for the potential. As shown in Fig. 1.1(a), the LJ potential can be divided

into the strong short-range repulsive component u0(r) and the slowly-varying long-

range attractive component u1(r):

u0(r) =


uLJ(r) + εLJ r < r0

0 r > r0

(1.2)

3



Figure 1.1: (a) The LJ potential and the its separation into the short and long-range
components based on WCA theory. (b) Comparison between the pair correlation
functions g(r) computed for the full LJ fluid and the WCA reference fluid at T ∗ =
0.65 and ρ∗ = 0.8442.

u1(r) =


−εLJ r < r0

uLJ(r) r > r0,

(1.3)

where r0 = 21/6σLJ is the distance where the LJ potential takes it minimum value.

The short-range component u0(r) captures the excluded volume repulsion between

LJ particles, and is referred to as the WCA potential. WCA theory suggests that

the vector sum of the long-range forces exerted on a given particle tends to cancel

at high density through the realization of force cancellation between oppositely

situated neighbors of the given molecule. Therefore, the average effect of u1(r)

is negligible in determining the structure of a dense uniform simple liquid, and

it merely contributes to a constant uniform background potential, allowing for a

mean filed theory treatment [13]. In contrast, the structure of the system can be

largely retrieved by the effect of molecular packing, which can be accounted for

by the excluded volume interaction captured by u0(r). This can be illustrated by

4



Figure 1.2: Comparison of the pair correlation functions for a large hard sphere
particle in the LJ fluid and for a large hard sphere particle in the WCA reference
fluid at T ∗ = 0.85 and ρ∗ = 0.72.

considering a WCA reference fluid composed of particles interacting only through

the pair potential u0(r) under the condition that it has the same density as the LJ

fluid. As shown in Fig. 1.1(b), the pair correlation function g(r) for the full LJ

fluid can be approximated by that of the WCA reference fluid to a good extent.

This assumption is utilized in the WCA perturbation theory for calculating the full

system’s thermodynamic properties such as the free energy [14].

When it comes to nonuniform fluids, the force cancellation argument no longer

holds. In an inhomogeneous or confined system, unbalanced forces emerge with the

introduction of large molecular interfaces. For instance, when a large hard sphere

solute particle is inserted into the LJ fluid, the long-range attractive forces on a LJ

particle near the solute-solvent interface do not cancel due to the large excluded vol-

ume created by the hard sphere solute. The interfacial structure associated with the

5



unbalanced long-range forces therefore cannot be captured by the WCA reference

fluid. This can be substantiated in Fig. 1.2 by the solute-solvent pair correlation

functions for the considered scenario. To achieve a physical description of the unbal-

anced forces in the inhomogeneous fluid systems, generalizations of the WCA ideas

are required. Local molecular field theory provides such a framework for treating

nonuniform fluids by formulating corrections to the WCA reference system.

1.2 Review of Local Molecular Field Theory

The structural and thermodynamic effects of the unbalanced forces emerging

in nonuniform fluids are accounted for in LMF theory by defining a mapping from

the full system described by a pair potential w(r) and an external potential ϕ(r) to

a mimic system described by the short-range component of w(r) and a renormalized

effective potential ϕR(r):

Full


w(r)

ϕ(r)

 LMF−−−→ Mimic


u0(r)

ϕR(r)

 ,

where w(r) = u0(r) + u1(r) follows the WCA idea of separating short and long-

range components of the intermolecular pair interaction. The short-range u0(r)

needs to be chosen on the appropriate molecular length scale such that it accurately

reproduces the local correlations in the full system by capturing the strong forces

between nearest-neighbor molecules. The average effects of the slowly-varying long-

range interactions are taken into account by the renormalized potential ϕR(r) as a

correction to the short-range reference in nonuniform systems.

6



The renormalized potential can be derived with a statistical mechanical treat-

ment utilizing the first equation of Yvon-Born-Green (YBG) hierarchy, which relates

the single-particle density ρ(r) to the pair density ρ(2)(r, r′). For the full system, it

reads

kBT∇ ln ρ(r) = −∇ϕ(r)−
∫

dr′ρ(r|r′)∇w(|r− r′|), (1.4)

where kB is the Boltzmann constant, and ρ(r|r′) is the conditional singlet density

defined by

ρ(r|r′) ≡ ρ(2)(r, r′)

ρ(r)
. (1.5)

It describes the density at r given that a particle is at r′. Similarly, for the mimic

system (denoted with the subscript R), the YBG equation reads

kBT∇ ln ρR(r) = −∇ϕR(r)−
∫

dr′ρR(r|r′)∇u0(|r− r′|). (1.6)

It is expected in LMF theory that, by appropriately choosing u0(r) and ϕR(r), the

nonuniform singlet density of the full system can be recovered by that of the mimic

system,

ρ(r; [ϕ]) = ρR(r; [ϕR]), (1.7)

where the functional dependence of the density on ϕ(r) or ϕR(r) is explicitly included

to facilitate subsequent analysis. With the assumption of Eq. (1.7), one can subtract

Eq. (1.4) from Eq. (1.6) to obtain

−∇ϕR(r) = −∇ϕ(r)−
∫

dr′ρR(r
′; [ϕR])∇u1(|r− r′|)

−
∫

dr′[ρ(r′|r; [ϕ])− ρR(r
′|r; [ϕR])]∇u0(|r− r′|)

−
∫

dr′[ρ(r′|r; [ϕ])− ρ(r′; [ϕ])]∇u1(|r− r′|).

(1.8)

7



A proper separation of short and long-range interactions can allow for further

simplification of Eq. (1.8). In the third term (line 2) of Eq. (1.8), as the gradient

of the short-range u0(r) exactly vanishes at large distances, this term only depends

on the difference between the short-range parts of conditional singlet densities for

the full and mimic systems. A proper choice of u0(r) guarantees that the two

systems encompasses the same strong short-range forces, which in turn control the

short-range part of the conditional densities. Consequently, the third term can

be neglected to a good approximation. In the fourth term (line 3) of Eq. (1.8),

provided that u0(r) captures the essential short-range strong forces, the gradient of

u1(r) is non-vanishing only at large distance, where the conditional singlet density

is approximately equal to the single particle density due to that systems away from

the critical point only possess short-range correlations. Therefore, the forth term

can be neglected as well. With the above approximations, the remaining terms in

Eq. (1.8) can be integrated to yield the LMF Equation:

ϕR(r) = ϕ(r) +

∫
dr′ρR(r

′; [ϕR])u1(|r− r′|) + C, (1.9)

where C is an integration constant. Eq. (1.9) is a self-consistent equation that can

be solved exactly by performing self-consistent loops between ϕR(r) and ρR(r). It is

derived with controlled and accurate approximations, and should not be considered

as a mean-field approach.

In the above context, LMF equation is derived for a single component system.

It can be readily applied to systems with charged sites as well [1]. In particular, the

Coulomb potential arising from charged particles decays at large distances with a

8



Figure 1.3: Separation of the Coulomb potential into the strong rapidly-varying
short-range component v0(r) and the slowly-varying long-range component v1(r),
with the smoothing length σ = 4.5 Å.

power smaller than the spatial dimension. Therefore any spherical truncation of it

leads to a diverging tail correction to the potential energy [16]. For the Coulomb

potential, the separation of short and long-range components can be performed as

v(r) ≡ 1

r
=

erfc (r/σ)

r
+

erf (r/σ)

r
≡ v0(r) + v1(r), (1.10)

where erf and erfc are the error function and complementary error function. The

long-range v1(r) arises from a unit Gaussian charge distribution with smoothing

length σ, defined as

ρG(r) =
1

π3/2σ3
exp

(
− r2

σ2

)
. (1.11)

As illustrated in Fig. 1.3, when the smoothing length σ is chosen to be on the

scale of characteristic nearest-neighbor distance of the system, v1(r) becomes suf-

ficiently slowly varying and approaches the Coulomb tail at distances greater than

the smoothing length. The strong short-range v0(r) is the screened potential due to

9



a point charge superposed with a neutralizing Gaussian charge distribution. It cap-

tures the forces from the Coulomb potential at distances less than σ while rapidly

vanishes at distances greater than σ.

Throughout this thesis, the LMF idea of separating short and long-rang in-

teractions will be utilized for constructing molecular models and analyzing the be-

haviors of nonuniform charged and polar aqueous systems in the relevant dynamic

processes.

1.3 Chapter Overview of Thesis

For the nonequilibrium process of solvation dynamics of ionic or dipolar ex-

citation in a solute molecule in water, the linear response relation between the

normalized dynamic Stokes shift and the normalized time correlation function of-

ten beaks down potentially due to the formation of hydrogen bonds between the

solute and water molecules, or due to the restructuring of water near the solute-

solvent interface. These effects can be attributed to the short-range solute-solvent

interactions. Therefore, the LMF idea of treating short and long-range interactions

separately can facilitate the examination of the origins of linear and nonlinear ef-

fects in solvation dynamics. In Chapter 2, an approach is formulated to relate the

dynamic linear response in ionic solvation dynamics to the long-range component

of solute-solvent Coulomb interactions in the aqueous environment. This approach

relies on the presumption that the slowly-varying electrostatic perturbation to the

bulk water induces minimal changes to the phase space distribution of the system.

10



Hence, for a Gaussian-smoothed ion excitation, the normalized dynamic Stokes shift

can be approximated by the normalized time correlation function obtained from

spontaneous fluctuations in equilibrium. Then it is further shown that the inclu-

sion of strong short-range solute-solvent interactions typically results in additional

nonlinear effects in solvation dynamics.

In Chapter 3, motivated by the conclusion from the previous chapter that

the linear effects in solvation dynamics are essentially associated with the slowly-

varying long-range electrostatic interactions, the connection between the intrinsic

dielectric dispersion of water and the solvation models idealized from LMF theory

is considered. The time scales in solvation dynamics are described differently for

ionic and dipolar solutes by the analytically theory that employs traditional hard

sphere solvation models. However, this time scale difference can be attributed to

the improper application of macroscopic electrostatic approaches to an atomically

sharp boundary. Using the Gaussian-smoothed charge and dipole models, it is shown

that the dynamic solvation response arising from the long-range electrostatic solute-

solvent interactions does not depend on the multipolar properties of the solute charge

distributions. This feature allows us to describe the time scales of solvation dynamics

in the linear response regime and to determine the dielectric dispersion intrinsic to

the solvent water. This chapter is based heavily on Renjie Zhao, Richard C Remsing,

and John D Weeks, Journal of Statistical Physics, 180:721–738, 2020 [17].

Chapter 4 and Chapter 5 concern a different dynamical process, crystal nucle-

ation, which occurs on a much greater time scale than solvation dynamics. Computer

simulations can provide us with the atomic scale information that are inaccessible

11



to experimental measurements in the course of nucleation. However, computational

studies of crystal nucleation in aqueous solutions are typically hindered by the ex-

traordinarily long time scale of the process relative to the conventional MD simu-

lations. This difficulty arises from both the computational overhead encountered

in treating Coulomb interactions in the aqueous environment and the rare event

nature of crystal nucleation. Therefore, LMF theory is introduced in this context to

reduce the model complexity for the long-range electrostatics, and it can be readily

combined with the advanced sampling methods that are widely used to accelerate

the sampling of rare events. In Chapter 4, the complex dynamic process of urea nu-

cleation in water is simulated by utilizing a LMF-based molecular model combined

with the advanced sampling method of well-tempered metadynamics [18]. Specif-

ically, LMF theory is firstly used to renormalize the effective long-range Coulomb

interactions between solute atomic sites. Secondly, well-tempered metadynamics

is employed to speed up rare events governed by short-range interactions. The

nucleation processes and the associated free energy landscape are successfully re-

produced through this combined approach. By comparing to the results simulated

with a short-range reference model, it can be observed that the contribution of the

long-range interactions to the stability of the the primary crystal state of urea in

water at ambient conditions.

In Chapter 5, the mechanism of urea nucleation in water is investigated in

detail. The one-step pathway described by classical nucleation theory is often vio-

lated for nucleation from aqueous solutions due to the existence of the intermediate

phase between the liquid phase and the crystal phase. Understanding the nucle-

12



ation mechanism for such cases requires an accurate description of the reaction

coordinates. To this end, the simulated trajectories of urea nucleation in water

are analyzed by an artificial intelligence-based approach that allows for accurately

approximating the reaction coordinates and separating the time scales of different

stages in the nucleation process. Through this analysis, the different roles of short

and long-range interactions in the nucleation process under consideration are dis-

cerned. It is revealed that the nucleation mechanism can largely be captured by

the short-range interactions, while the long-range interactions are indispensable in

accounting for metastability. The direct evidence for the two-step nucleation mech-

anism is also highlighted in the results. Chapter 4 and Chapter 5 are based heavily

on Renjie Zhao, Ziyue Zou, John D Weeks, and Pratyush Tiwary, arXiv preprint

arXiv:2307.13189, 2023 [19].

There are other types of physical processes in the aqueous environment in

addition to the above two topics discussed in this thesis. However, short and long-

range interactions do not always have distinct and non-trivial effects on a dynamical

process. For instance, in a component of my doctoral research work, it is found that

the long-range component of Coulomb interactions barely have any influence on

the pairing process of monovalent ions [20]. LMF-based models characterizing the

effects of long-range interactions are not utilized for this system. Therefore, this

paper is not closely connected to the theme of the thesis and is not included in the

thesis.

13



Chapter 2: Solvation Dynamics in Water I: Origins of Linear and

Nonlinear Dynamic Responses from the Local Molecular

Field Theory Perspective

2.1 Introduction

Solvation dynamics deals with the nonequilibrium response of the solvent en-

vironment to the redistribution of charge densities in a solute. This process, phys-

ically realized by electronic excitations or electron transfer reactions, can influence

the efficiency of relevant chemical and biochemical reactions [21–26]. In the typical

scenario of solvation dynamics where a point charge or dipole moment is instanta-

neously created inside an initially neutral solute core, the electrostatic interaction

energy E(t) between the excited charge density and the solvent molecules serves as

the natural observable during the equilibration of the system. In experiments, the

characteristics of energy relaxation is obtained by measuring the time-dependent

shift of the fluorescence emission maximum (Stokes shift) [27–33], and the results

are expressed through the normalized dynamic Stokes shift function,

S(t) =
∆E(t)−∆E(∞)

∆E(0)−∆E(∞)
, (2.1)

14



which is termed as normalized nonequilibrium response function to accommodate the

theoretical context here. In Eq. (2.1), ∆E(t) denotes the nonequilibrium ensemble

average of the change in the electrostatic interaction energy upon the redistribution

of solute charge density.

In computational studies, there have been efforts [34–36] to approximate the

normalized nonequilibrium response function S(t) by the normalized time correla-

tion function (TCF) in equilibrium

C(t) =
⟨δE(t)δE(0)⟩
⟨δE(0)δE(0)⟩

, (2.2)

where the brackets ⟨⟩ denotes the ensemble average of the system, and δE(t) denotes

the fluctuation of the electrostatic interaction energy from its equilibrium ensemble

average at time t. The basis of this approximation can be traced back to Onsager’s

regression hypothesis, which suggests “the average regression of fluctuations will

obey the same laws as the corresponding macroscopic irreversible processes” [37],

i.e. the relaxation of spontaneous microscopic fluctuations in a system evolves with

time in the same manner as the nonequilibrium response of the system to small

perturbations. The detailed validation for the approximation of S(t) by C(t) re-

lies on the applicability of linear response theory to the specific solvation dynamics

problems under consideration, and it has been observed that linear response theory

breaks down in some instances of solvation dynamics problems [34, 38–40]. The

traditional approach to linear response theory [35, 41] requires the perturbation to

the system H′ being sufficiently small relative to kBT , such that the Boltzmann fac-

tor of the perturbation can be approximated by its first-order expansion, 1 + βH′.

15



However, there are circumstances in solvation dynamics where this assumption can-

not be consistently satisfied. The electrostatic interaction energy due to the newly

created ion or dipole in the solvent can often exceed kBT and, in some cases, even

surpass it by two orders of magnitude. To overcome this limitation, a more general

condition was proposed for the linear response approximation, where it assumes

that the equilibrium fluctuation of the solvation energy δE(t) in the excited state is

a Gaussian random variable [42, 43]. Nonetheless, this condition only concerns the

statistical characteristics of the linear response, and the molecular details related to

linear and nonlinear responses in solvation dynamics are yet fully understood.

In this chapter, the molecular origins of the linear and nonlinear responses in

solvation dynamics are analyzed and computationally examined from the perspec-

tive of LMF theory. For the solvation of a ionic solute in equilibrium, the solvation

free energy is composed of the electrostatic and non-electrostatic components. The

non-electrostatic component of the solvation free energy arises from the repulsive

excluded volume interactions and the van der Waals interactions between the solute

core and the solvent molecules even before the existence of solute charges. In the

aqueous environment, the introduction of the solute core induces a molecular scale

interface that perturbs the local hydrogen-bond network of the bulk water. As a re-

sult, non-zero solvent charge densities are formed around the solute-solvent interface

and it further contributes to the electrostatic response when the solute charges are

presented [44]. This contribution has been shown to be nonlinear and asymmetric

with respect to the sign of the solute charge in the charging component of the ionic

solvation free energy [45]. It is therefore natural to speculate that in situations where

16



the bulk water structure is substantially perturbed by the solute core, the nonlinear

and asymmetric effects resulting from the solute core in equilibrium would manifest

in solvation dynamics as well.

Apart from the effects induced by the solute core, ion solvation also encom-

passes the response of the solvent to the Coulomb potential arising from the solute

charge. LMF theory separates the Coulomb potential into the strong short-range

component v0(r), which is rapidly varying over molecular length scales, and the

slowly-varying long-range component v1(r) as

v(r) ≡ 1

r
=

erfc (r/σ)

r
+

erf (r/σ)

r
≡ v0(r) + v1(r), (2.3)

where erf and erfc are the error function and complementary error function, respec-

tively. While a portion of the short-range component, particularly the diverging

part as r approaches zero, is screened by the solute core, solvent molecules are still

subjected to the influence of this strong short-range electrostatic interaction. The

rapidly-varying nature of the short-range electrostatic interaction typically leads

to induced solvent charge densities that are nonlinear with respect to the solute

potential [45].

To isolate the anticipated linear response in ion solvation, we consider the

long-range electrostatic interaction Qv1(r), which can be equivalently interpreted as

arising from the physical model of a Gaussian charge distribution. Depending on

magnitude of charge Q, one can choose the smoothing length l of the Gaussian on the

scale of nearest-neighbor correlations of the system to ensure Qv1(r) is sufficiently

slowly varying, such that the solvent response to the Gaussian charge distribution

17



obeys linear response. Inserting the Gaussian charge distribution into the solvent

without invoking the solute core would not pose any conceptual difficulty because its

potential remains finite at r = 0. In this manner, one is able to calculate the linear

and symmetric component (and in general the dominant component) of the charging

free energy in ion solvation [45]. In the context of solvation dynamics, it is also of

interest to investigate whether the response to a Gaussian charge distribution follows

linear response theory. The condition of the perturbation potential being slowly

varying over space, as satisfied by a Gaussian charge distribution with a suitably

chosen LMF length scale l, is equivalent to keeping only small k components in

the Fourier space. It is a physically more general and intuitive condition for linear

response, in contrast to the traditional approach that relies on the perturbation

Hamiltonian being small on the scale of kBT .

2.2 Linear Response to Gaussian Charge Distributions in Equilib-

rium

In this section, it is shown that in equilibrium systems, charge densities which

generates slowly-varying electrostatic potential induce linear responses in the sol-

vent. Consider inserting a charge distribution ρQ(r) in a uniform and isotropic liquid

system. The potential arising from ρQ(r) is

vQ(r) =

∫
ρQ(r′)

|r− r′|
dr′. (2.4)

18



The interaction energy between the charge distribution ρQ(r) and the solvent for a

particular microscopic molecular configuration R is

E(R) =

∫
ρq(r,R)vQ(r)dr, (2.5)

where

ρq(r,R) =

NC∑
i=1

qiδ(r− ri(R)) (2.6)

is the solvent charge density in configuration R, which is composed of NC sites

each with charge qi located at position ri(R). The average solvent charge density

ρq(r) can be obtained by taking the ensemble average in the presence of the inserted

charge distribution ρQ(r),

ρq(r) = ⟨ρq(r,R)⟩Q. (2.7)

By treating vQ(r) as an weak external perturbation to the liquid system, func-

tional derivatives of the grand potential Ω [46] yield

δ[βΩ]

δ[−βvQ(r)]
= ρq(r) (2.8)

and

δ2[βΩ]

δ[−βvQ(r)]δ[−βvQ(r′)]
=

δρq(r)

δ[−βvQ(r′)]
= χqq(r, r′), (2.9)

where χqq(r, r′) ≡ ⟨δρq(r,R)δρq(r′,R)⟩0 is the charge-charge correlation function

evaluated as the ensemble average without the presence of the external perturbation

potential, with the fluctuation in the solvent charge density denoted by δρq(r,R).

Based on Eq. (2.9), we can approximate the induced solvent charge density by a

linear functional of the external potential,

ρq(r) ≈ −β

∫
dr′χqq(r, r′)vQ(r′), (2.10)

19



where the average solvent charge density is taken to be zero in the absence of vQ(r)

because of neutrality.

Using the approach of coupling parameter method [41, 45], we assign a linear

coupling parameter λ to the external potential, i.e. vQλ (r) ≡ λvQ(r). Then the Gibbs

free energy of a system coupled with λ follows

Gλ ∝ −kBT ln

∫
exp

[
− β

(
H0 + Eλ(R)

)]
dR, (2.11)

where

Eλ(R) =

∫
ρq(r,R)vQλ (r)dr, (2.12)

and H0 is the Hamiltonian of the unperturbed system. Differentiating the Gibbs

free energy yields

dGλ

dλ
=
〈dEλ(R)

dλ

〉
λ

=

∫
ρqλ(r)v

Q(r)dr,

(2.13)

where ρqλ(r) =
〈
ρq(r,R)

〉
λ
. By integrating Eq. (2.13), the free energy due to the

external potential can be obtained as

∆G =

∫ 1

0

dλ

∫
drρqλ(r)v

Q(r). (2.14)

Presuming that linear response approximation is valid, we insert Eq. (2.10) into Eq.

(2.14), yielding

∆G = −β

∫ 1

0

dλ

∫
drvQ(r)

∫
dr′χqq(r, r′)vQλ (r

′)

= −β

∫ 1

0

λdλ

∫∫
drdr′vQ(r)χqq(r, r′)vQ(r′)

= −β

2

∫∫
drdr′vQ(r)

〈
δρq(r,R)δρq(r′,R)

〉
0
vQ(r′)

= −β

2

〈
E2(R)

〉
0
.

(2.15)

20



This result can be alternatively derived by assuming the fluctuation in the interac-

tion energy follows Gaussian statistics. In the framework of thermodynamic pertur-

bation theory [47–49], the exact expression of the free energy due to the external

potential

∆G = −kBT ln
〈
exp

[
− βE(R)

]〉
0

(2.16)

takes form of a moment-generating function [50]. Therefore, in the equilibrium

ensemble without the presence of the external potential, if the interaction energy

E(R) is a Gaussian random variable, then the cumulant expansion of Eq. (2.16)

only includes the terms of the first two orders, i.e.

∆G =
〈
E(R)

〉
0
− β

2

〈(
E(R)− ⟨E(R)⟩0

)2〉
0
, (2.17)

and all higher order terms are necessarily zero. Since the bulk liquid system has zero

average charge density throughout space, we have ⟨E(R)⟩0 = 0. As a consequence,

Eq. (2.17) becomes identical to Eq. (2.15).

Having established the equivalence between Gaussian statistics and linear re-

sponse theory in equilibrium, we are now poised to investigate the suitability of

linear response theory in the context of Gaussian charge solvation. Fig. 2.1 depicts

the probability distribution of E(R) in equilibrium as if the system were subjected

to the influence of a Gaussian charge distribution with smoothing length l = 6 Å.

The statistics can be accurately described by a Gaussian random variable. We can

therefore presume the response of the bulk water system to a Gaussian charge dis-

tribution with an appropriate smoothing length follows linear response theory in

equilibrium. Indeed, in Ref. [45] it has been shown that the solvation free energies

21



Figure 2.1: Probability distribution (blue dots) of the interaction energy E(R)
between a Gaussian charge distribution with smoothing length l = 6 Å and the
bulk water modelled by SPC/E water, drawn from MD simulations. The simulation
results are fitted by a Gaussian process (black solid line).

arising from Gaussian charge distributions with adequately large smoothing lengths

l can be predicted by an analytical theory based on linear response.

2.3 Linear Response Theory for Solvation Dynamics

Linear response theory for nonequilibrium processes concerns the time evo-

lution of a system under the influence of a time-dependent perturbation weakly

coupled to the system. In contrast to the situations in equilibrium where the linear

response of the system is related to a category of density-density correlation func-

tions, in dynamic processes, the linear response of a system can be characterized

through time correlation functions that reflect the system’s spontaneous fluctua-

tions, leading to the well-known linear response relation [41].

Here we focus on the linear response approach to solvation dynamics. Unlike

22



the traditional approach in which assumptions about the magnitude of the pertur-

bation are made [35,41], or the more recent theory which requires the fluctuations in

the excited sate obeying Gaussian statistics [42], we follow Kubo’s approach [46,51]

to formally derive the linear response in solvation dynamics. We argue that it is the

property of the external potential that determines the system’s behavior in dynamic

responses, and the Gaussian charge distribution is the physical entity associated with

the linear component in solvation dynamics of ion.

By treating the excitation of the solute charge density as an external perturba-

tion potential changing with time vQ(r, t), the Hamiltonian of the system is denoted

as

H = H0 +H′(t), (2.18)

where H0 is the Hamiltonian of the unperturbed system, and H′(t) is the time-

dependent solute-solvent electrostatic interaction energy

H′(t) =

∫
ρq(r,R; t)vQ(r, t)dr

≡ ρq(r,R; t) · vQ(r, t),
(2.19)

where for the convenience of expression we denote the integral of a function with

vQ(r, t) over space by the inner product “ ·vQ(r, t)”. The time argument in ρq(r,R; t)

is extracted from the implicit time dependence in R, and is omitted for time t.

Consider the phase space distribution function of the system f(t) ≡ f(rN ,pN ; t),

23



the evolution of which is governed by the Liouville equation as

∂f(t)

∂t
= −iLf(t)

= {H0 +H′, f(t)}

= −iL0f(t) + {ρq(r,R), f(t)} · vQ(r, t),

(2.20)

where L is the Liouville operator defined as

L ≡ i{H, }, (2.21)

and {, } represents the Poisson bracket defined as

{A,B} ≡
∑
i

(
∂A

∂ri

∂B

∂pi

− ∂A

∂pi

∂B

∂ri

)
, (2.22)

where ri denotes the canonical coordinates of the system and should not be confused

with the argument r in ρq(r,R) and vQ(r, t). Assuming that at time t → −∞, the

system is in a thermal equilibrium state, i.e. f(−∞) = f0, we can rewrite f(t) as

f(t) = f0 +∆f(t), (2.23)

and by definition we have

∂f0
∂t

= −iL0f0. (2.24)

Inserting Eq. (2.23) and Eq. (2.24) into Eq. (2.20) and expanding the equation to

the first order in vQ(r, t), we obtain

∂∆f(t)

∂t
= −iL0∆f(t) + {ρq(r,R), f0} · vQ(r, t). (2.25)

The above linearization leading to Eq. (2.25) is based on the presumption that

∆f(t), the change in the phase space distribution due to the perturbation, has

24



small partial derivatives with respect to the canonical coordinates. The solution to

Eq. (2.25) is

∆f(t) =

∫ t

−∞
exp [−i(t− s)L0]{ρq(r,R), f0} · vQ(r, t)ds. (2.26)

Since the thermal state f0 ∝ exp (−βH0), we have{
ρq(r,R), f0

}
=
∑
i

(
∂ρq(r,R)

∂ri

∂f0
∂pi

− ∂ρq(r,R)

∂pi

∂f0
∂ri

)

= −β
∑
i

(
∂ρq(r,R)

∂ri

∂H0

∂pi

− ∂ρq(r,R)

∂pi

∂H0

∂ri

)
f0

= −β(iL0ρ
q(r,R))f0

= −βρ̇q(r,R)f0.

(2.27)

With Eq. (2.26) and Eq. (2.27), we then examine the evolution of the nonequi-

librium ensemble average of ρq(r,R) upon the influence of the external potential,

which reads

ρq(r; t) =

∫∫
ρq(r,R)f(t)drNdpN

=

∫∫
ρq(r,R)∆f(t)drNdpN + ρq0(r)

= −β

∫ t

−∞
ds vQ(r′, s) ·

∫∫
f0ρ

q(r,R) exp [−i(t− s)L0]ρ̇q(r
′,R)drNdpN

= −β

∫ t

−∞
ds vQ(r′, s) ·

∫∫
f0ρ̇q(r

′,R) exp [i(t− s)L0]ρ
q(r,R)drNdpN

= −β

∫ t

−∞

〈
ρq(r,R; t− s)ρ̇q(r′,R; 0)

〉
0
· vQ(r′, s)ds

= β

∫ t

−∞

〈
ρ̇q(r,R; t− s)ρq(r′,R; 0)

〉
0
· vQ(r′, s)ds,

(2.28)

where we have dropped ρq0(r) in the second line of the equation because of neutrality.

Eq. (2.28) describes the linear relation between the time-dependent response in the

25



solvent charge density and the external perturbation potential. In the situations of

ion excitation, the change in the solute charge density takes place instantaneously.

Therefore, we prescribe that the time dependence of the external potential takes the

form

vQ(r, t) = vQ(r)Θ(t), (2.29)

where Θ(t) is the Heaviside step function. Inserting Eq. (2.29) into the result of

Eq. (2.28), we have

ρq(r; t) = β

∫ t

0

〈
ρ̇q(r,R; t− s)ρq(r′,R; 0)

〉
0
· vQ(r′)ds (2.30)

Make the change of variable of integration from s to t− s and we obtain

ρq(r; t) = −β

∫ t

0

〈dρq(r,R; t− s)

d(t− s)
ρq(r′,R; 0)

〉
0
· vQ(r′)d(t− s)

= β

(〈
ρq(r,R; t)ρq(r′,R; 0)

〉
0
−
〈
ρq(r,R; 0)ρq(r′,R; 0)

〉
0

)
· vQ(r′).

(2.31)

Then taking the operation vQ(r)· on both sides of Eq. (2.31) yields

E(t) = vQ(r) · ρq(r; t)

= βvQ(r) ·

(〈
ρq(r,R; t)ρq(r′,R; 0)

〉
0
−
〈
ρq(r,R; 0)ρq(r′,R; 0)

〉
0

)
· vQ(r′)

= β

(〈
vQ(r) · ρq(r,R; t)ρq(r′,R; 0) · vQ(r′)

〉
0

−
〈
vQ(r) · ρq(r,R; 0)ρq(r′,R; 0) · vQ(r′)

〉
0

)

= β
(
⟨E(t)E(0)⟩0 − ⟨E2(0)⟩0

)
,

(2.32)

where the electrostatic interaction energy E(t) is given by the perturbation Hamilto-

nian Eq. (2.19). Eq. (2.32) is a special case of the fluctuation-dissipation theorem,

26



through which we can relate the normalized nonequilibrium response function S(t)

to the normalized TCF:

S(t) =
∆E(t)−∆E(∞)

∆E(0)−∆E(∞)

=
β
(
⟨E(t)E(0)⟩0 − ⟨E2(0)⟩0

)
− β

(
⟨E(∞)E(0)⟩0 − ⟨E2(0)⟩0

)
β
(
⟨E(0)E(0)⟩0 − ⟨E2(0)⟩0

)
− β

(
⟨E(∞)E(0)⟩0 − ⟨E2(0)⟩0

)
=

⟨E(t)E(0)⟩0 − ⟨E(∞)E(0)⟩0
⟨E(0)E(0)⟩0 − ⟨E(∞)E(0)⟩0

.

(2.33)

Note that in the limit of t → ∞, E(t) becomes uncorrelated to E(0). Therefore, we

have

lim
t→∞

⟨E(t)E(0)⟩0 = ⟨E⟩20, (2.34)

where ⟨E⟩20 is the square of the equilibrium interaction energy in the ensemble

without the presence of the perturbation potential. Based on this observation we

further have

S(t) =
⟨E(t)E(0)⟩0 − ⟨E⟩20
⟨E(0)E(0)⟩0 − ⟨E⟩20

=
⟨δE(t)δE(0)⟩0
⟨δE(0)δE(0)⟩0

= C(t).

(2.35)

2.4 Molecular Origins of Nonlinear Effects

The formal derivation of the equivalence between S(t) and C(t) relies on the

assumption that the potential due to the excited charge density induces linear re-

sponses in the systems, such that the linearized equation Eq. (2.25) holds true. This

assumption can be fulfilled when dealing with a Gaussian charge distribution due to

the slowly-varying nature of its electrostatic potential. When considering Gaussian

27



charge distributions with relatively modest magnitude of the total charge (Q < 4e0),

choosing the Gaussian width on the scale of typical nearest neighbor distances in

the solvent (l > 3 Å) can ensure the validity of linear response theory under equi-

librium conditions [45, 52]. It is therefore natural to anticipate the excitation of a

Gaussian charge distribution with the width on a comparable scale would also lead

to solvation dynamics that follows linear response in nonequilibrium scenarios.

Figure 2.2: (a) Comparison between S(t) (the normalized nonequilibrium response
function) and C(t) (the normalized TCF) for the excitation of a Gaussian charge
distribution with Q = 3e0 and l = 6 Å in SPC/E water. (b) Comparison between
S(t)’s for the excitation of Gaussian charge distributions with Q = ±3e0 and l =
6 Å in SPC/E water.

Here we monitor the solvation dynamics for an instantaneously excited Gaus-

sian charge distribution with Q = 3e0 and l = 6 Å in a bulk water system computed

with MD simulations. The Gaussian charge distribution is realized as an external

electrostatic potential, and the solvent is modelled by SPC/E water. Fig. 2.2(a)

shows that the S(t) for the Gaussian charge distribution can be well approximated

by the normalized TCF, with merely negligible differences. Furthermore, when con-

sidering two Gaussian charge distributions with charges of equal magnitudes but

28



differing in sign, their dynamics yield indistinguishable S(t) profiles, as illustrate in

2.2(b). This suggests that the solvation dynamics for Gaussian charge distributions

within the relevant parameter range not only follows linear response, but also ex-

hibits symmetry with respect to the sign of charge, a noteworthy feature previously

demonstrated in the equilibrium contexts [45].

However, in realistic models of ionic solvation dynamics, additional molecular

details associated with nonlinear effects generally exist. In classical models of ion,

the solute charge is represented by a point charge embedded in a solute core. The

Coulomb interaction arising from the point charge includes the short-range compo-

nent as well as the long-range component which is accounted for by the Gaussian

charge distribution with the properly chosen smoothing parameter. While the solute

core does screen a fraction of the short-range component of the Coulomb interac-

tion, the remaining strong short-range interaction still exerts its influence on solvent

molecules to different extents depending on the size of the solute core. To exam-

ine the effect of the short-range electrostatic interaction in solvation dynamics, we

consider a Gaussian charge distribution with the same total charge Q = 3e0 and a

smaller smoothing length l = 2 Å. This is amount to including some larger k com-

ponents in the external electrostatic potential, as compared to the previous case.

In this phase, even though we are not aiming to take into account the impact of

the solute core, the Gaussian width of l = 2 Å is chosen such that, in the spatial

range r ≥ 3 Å, the potential generated by the Gaussian charge closely approxi-

mates that of a point charge situated within a solute core with an effective radius

of Rc = 3 Å (as illustrated by Fig. 2.3(a)). Since the Gaussian charge behaves con-

29



sistently throughout space, the situation under consideration guarantees that the

Gaussian charge density exhibits a spatial rate of variation in electrostatic potential

akin to that observed for a realistic solute model with a solute core.

Figure 2.3: (a) Comparison between the Coulomb potential (black solid line) and
the potential v1(r) (red dotted line) arising from the Gaussian charge distribution
with l = 2 Å. Since the smoothing length l = 2 Å is not large enough for the
potential v1(r) to be slowly varying, the plotted v1(r) contains short-range electro-
static interactions on the relevant molecular length scale. The violet auxiliary line
at r = 3 Å indicates a potentially existing solute core with an effective radius of
3 Å. In the range r ≥ 3 Å, v1(r) closely approximates the Coulomb potential. (b)
Comparison between S(t) and C(t) for the excitation of a Gaussian charge distri-
bution with Q = 3e0 and l = 2 Å in SPC/E water.

Fig. 2.3(b) shows the solvation dynamics for the excitation of such a Gaussian

Charge distribution with Q = 3e0 and a relatively small smoothing length of l =

2 Å. Evident inconsistencies are observed between the profiles of S(t) and C(t),

especially during the oscillatory period of energy relaxation, in notable contrast to

the case shown in Fig. 2.2(a). It is the utilization of the smaller smoothing length

l that leads to the deviation of the electrostatic perturbation potential from being

adequately slowly varying, the condition we presume necessary for linear response

theory to hold in the current context. Consequently, we can ascribe the above

30



observed nonlinear effects in ionic solvation dynamics to the response elicited by the

short-range component of the electrostatic potential arising from the excited solute

charge.

Naturally a neutral solute core always preexists the excitation of any solute

charge densities. This excluded volume core can strongly perturb the uniformity of

the bulk water system, resulting in significant local structure rearrangements and

non-zero charge densities ρq0(r) near the solute-solvent boundary [45,52], which may

in turn induce nonlinear responses upon the excitation of the solute charge. Here,

in order to resemble the realistic situation of ionic solvation dynamics, we prepare a

neutral solute core in equilibrium with the water system prior to exciting the solute

charge density. The solute core is modelled by a excluded volume hard sphere with

a radius Rc = 3 Å for the sake of simplicity. However, the incorporation of van

der Waals interactions would not qualitatively change the behaviors of the system.

Then we instantaneously create a point charge with Q = 3e0 centred with the hard

sphere.

Fig. 2.4(a) displaces the solvation dynamics for the process described above.

With the introduction of the influence from the solute core, the distinctions between

the profiles of S(t) and C(t) have become more pronounced and complicated. The

differences between S(t) and C(t) in the oscillatory period of energy relaxation are

more remarkable than in the situation of bare Gaussian charge excitation shown in

Fig. 2.3(b). Moreover, S(t) and C(t) exhibit quite different relaxation rates after

the the oscillatory period. The asymptotic relaxation rate in solvation dynamics is

related to the dielectric properties in the high frequency limit and in the equilibrium

31



Figure 2.4: (a) Comparison between S(t) and C(t) for the excitation of a point
charge with Q = 3e0 embedded in a hard solute core of Rc = 3 Å in SPC/E water.
(b) Comparison between S(t)’s for the excitation of point charges with Q = ±3e0
embedded in a hard solute core of Rc = 3 Å in SPC/E water.

limit [17], and this topic will be detailed in the next chapter. It is therefore indi-

cated that the existence of the solute core has impacted on the dielectric response

of the system. Such an effect is reasonably anticipated, considering that the solute

core disrupts local charge neutrality even prior to the excitation of the solute charge

distribution. Fig. 2.4(b) compares the actual dynamics, represented by S(t), be-

tween the solvation processes upon exciting the point charges with Q = ±3e0 in the

presence of the solute core. The inconsistency between the two profiles indicates an

asymmetric dynamic effect with respect to the sign of charge. The above observa-

tions suggests that the preexisting solute core introduces additional highly nonlinear

effects to solvation dynamics, as compared to the case of exciting a Gaussian charge

in the uniform water system, and these nonlinear effects are also asymmetric with

respect to the sign of charge. It should be noted that other kinetic factors, such

as the initial momentum and angular momentum of the solute, may also affect the

behavior of the system in terms of nonlinearity [53,54], but they are not the focus of

32



this work. In the MD simulations, the position of the solute is fixed and its kinetic

energy is neglected.

2.5 Conclusion

In this chapter, both linear and nonlinear effects are investigated in the context

of ionic solvation dynamics in water, adopting the perspective of LMF theory. We

evaluate the applicability of linear response theory to solvation dynamics in various

scenarios by examining the proximity between the normalized nonequilibrium re-

sponse function S(t) and the normalized Time Correlation Function C(t), obtained

from MD simulations.

Facilitated by the model of the Gaussian charge distribution, we effectively

dissect the different molecular details within the solvation model that pertain to

linear and nonlinear effects. LMF theory separates intermolecular interactions into

short and long-range components. A Gaussian charge distribution with the smooth-

ing parameter appropriately chosen on the relevant molecular length scale generates

the electrostatic potential equivalent to the long-range component of the Coulomb

interaction arsing from a realistic charged solute. We review the linear response

theory for the solvation of Gaussian charge distributions in equilibrium, and note

that it is the slowly-varying nature of the long-range interaction that enables linear

response in this context. Then we extend this approach to the nonequilibrium sce-

narios, where the interaction energy of an aqueous system relaxes subsequent to the

excitation of the solute charge density.

33



Through simulating the excitation of Gaussian charge density in the bulk wa-

ter, it is revealed that the response to the long-range electrostatic potential arising

from the solute charge accurately follows linear response theory in solvation dy-

namics. However, the inclusion of the strong short-range electrostatic interaction

between the solute and the solvent in this process results in nonlinear dynamic re-

sponses, clearly evidenced by the breakdown of the linear response relation. More-

over, we account for the preexisting solute core, which can be classified as a potent

short-range intermolecular interaction within the LMF framework. This element

introduces additional non-trivial nonlinear and asymmetric effects into solvation

dynamics.

Having established the connection between the linear effects in solvation dy-

namics and the physical model of Gaussian charge distribution with slowly-varying

long-range electrostatic potential, it will be further investigated in the next chapter

how the Gaussian charge model reveals the intrinsic dielectric dispersion of water in

solvation dynamics.

34



Chapter 3: Solvation Dynamics in Water II: Intrinsic Dielectric Dis-

persion of Water1

3.1 Introduction

As solvation dynamics captures the temporal polarization response of a po-

lar solvent to changes in charge distribution, solvation models can serve as natu-

ral probes for the dielectric dispersion ε(ω), i.e the frequency-dependent dielectric

property, of the solvent. Classic models of ionic and dipolar solvation are based on

dielectric continuum theory, in which the polar solvent is treated as a uniform fea-

tureless dielectric medium with dielectric permittivity ε. The dielectric property of

the solvent can be probed by inserting an approximation to the solute charge to the

uniform dielectric solvent. The theoretical difficulty arises as an idealized point test

charge or dipole generates the singularity in the electrostatic potential in the limit

of r → 0, which when overlapping with the solvent charges, results in an infinite

interaction energy. To overcome this issue, classic approaches employ an excluded

volume solute core to shield its embedded point charge or dipole from the charges

in the solvent.

1This chapter and Appendix A are based heavily on ‘Renjie Zhao, Richard C Remsing, and

John D Weeks. Journal of Statistical Physics, 180:721–738, 2020’ [17].

35



The Born and Bell models [55, 56] are the most typical representatives of the

dielectric continuum theory-based solvation models described above. They are con-

structed by placing a point charge or point dipole in the centre of a spherical hard

solute core. Though the Born-like models are utilized to estimate the charging

solvation free energy to the correct order of magnitude, their absence in address-

ing solvent structuring at the boundary of the solute core can result in qualitative

inaccuracies related to the nonlinear and asymmetric solvent response to the phys-

ically realistic solute model [45]. Moreover, it is well known that the Born and Bell

models predict different relaxation behaviors (in particular the solvent longitudinal

relaxation time) of the systems depending on the order of the multipole embedded

in the solute core [57–59]. Therefore, the question arises as to whether these dif-

ferent dynamic behaviors stem from the multipolar properties of the solute charge

distributions.

Here, to investigate the dielectric responses of ionic and dipolar solvation,

we first develop a theory for the solvation of Gaussian dipoles in polar solvents

as the Gaussian counterpart to the Bell model. Similar to the case of Gaussian

charge distributions, by choosing the smoothing length to be on the order of molec-

ular correlations, the Gaussian dipole model captures the slowly-varying long-range

component of the electrostatic potential arising from an ideal dipole. Therefore,

its interaction energy with atomically-detailed solvent models can be described by

linear response theory to quantitative accuracy. The utilization of Gaussian charge

and Gaussian dipole allows us to to rule out the nonlinear effects caused by the

short-range solute-solvent interactions, which as we have shown in Chapter 2, can

36



altered the detailed relaxation rate in solvation dynamics.

The analytical descriptions developed for solvation dynamics of Gaussian charge

and dipole distributions in the framework of linear response theory are then com-

pared to the results of MD simulations. We demonstrate that as the Gaussian charge

and Gaussian dipole distributions account for the linear components in the solva-

tion responses of monopole and dipole, respectively, their consistent linear dynamic

responses do not depend on their multipolar properties. Therefore, in this man-

ner, we are able to isolate the dielectric response intrinsic to the dielectric property

of the solvent, without implicit assumptions of the nature of the solvent response

to solute cores. Finally, the dielectric dispersion of a widely-used water model is

determined by analyzing the simulation results of solvation dynamic responses of

Gaussian charge and dipole distributions in the framework of linear response theory.

3.2 Solvation of Dipolar Solutes in Equilibrium

3.2.1 Bell Model

The Bell model, categorized as a Born-like model for describing dipole solva-

tion, is constructed by placing a point dipole in the centre of an excluded volume

hard sphere with radius R, while the solvent is idealized as the featureless dielectric

continuum outside the hard sphere [56]. Similar to the Born model, the conventional

treatment for the Bell model assumes that the system responds linearly to the elec-

tric dipole moment. Moreover, the manner in which the interaction energy or the

charging free energy of the Bell model depends on the dielectric permittivity ε has

37



further consequences on solvation dynamics. In this section, we begin by revisiting

this model before dealing with its Gaussian counterpart.

Continuing with the coupling parameter method discussed in Section 2.2, the

charging free energy is calculated with a charging process in which the solute charge

distribution evolves as a partially charged point dipole that generates the electro-

static potential vµλ (r) = λvµ(r), where vµ(r) is the electrostatic potential of a point

dipole with a dipole moment along the z-axis of magnitude µ,

vµ (r) =
µ cosφ

r2
, (3.1)

and λ is a linear coupling parameter varying between 0 and 1. Analogous to the

coupling parameter integration of ionic solvation, Eq. (2.14), the integral for the

charging component of the solvation free energy of dipole takes the form

∆G (µ) =

∫ 1

0

dλ

∫
drρqλ(r)v

µ(r), (3.2)

where the Bell model, by definition, assumes zero preexisting solvent charge density

arising from hypothetical and/or structural boundaries, i.e. an infinite solvent is

considered [52].

The superposed electrostatic potential due to both the solute charge distribu-

tion and the induced solvent charge density can be calculated by solving Poisson’s

equation with the relevant boundary condition

vBell
λ (r, φ) =


3

2ε+ 1

λµ cosφ

r2
r > R

λµ cosφ

r2
− 2 (ε− 1)

2ε+ 1

λµr cosφ

R3
r < R ,

(3.3)

where ε is the dielectric permittivity of the solvent and we have assumed that the

dipole is aligned long the z-axis, µ = µẑ, without loss of generality. The induced

38



solvent charge density ρqλ(r) on the surface of the solute core is then given by

4πρqλ (r) = −∂vBell
λ

∂r

∣∣∣∣
r=R+ϵ

+
∂vBell

λ

∂r

∣∣∣∣
r=R−ϵ

(3.4)

which yields

ρqλ (r) = −3λµδ(r −R) cosφ

2πr3
ε− 1

2ε+ 1
. (3.5)

Provided that the induced charge density of the Bell model depends linearly on the

coupling parameter λ, we can integrate out λ and Eq. (3.2) reduces to

∆G (µ) =
1

2

∫
drvµ (r) ρq (r), (3.6)

where ρq (r) = (r) ρqλ=1 is the induced charge density of the fully coupled solute-

solvent system.

Inserting Eq. (3.5) and Eq. (3.1) into Eq. (3.6), we can obtain the charging

free energy of the Bell model

∆G (µ) = −3µ2

4π

ε− 1

2ε+ 1

∫ 2π

0

dθ

∫ π

0

dφ cos2 φ sinφ

∫ ∞

R

dr
δ (r −R)

r3

= − ε− 1

2ε+ 1

µ2

R3
.

(3.7)

This result can be equivalently derived from the reaction field approach [60], which

evaluates the work needed for bringing the solute charge distribution from infinitely

far away to its actual position with infinitesimally small steps under the influence

of the the polarization field (reaction filed) arsing from the induced solvent charge

density.

Similar to the solvent charge density of the Born model, the solvent charge den-

sity of the Bell model is localized on the surface of the solute core. The discontinuity

39



in this charge density in the radial direction results from applying macroscopic elec-

trostatic treatments to an atomically sharp boundary. Moreover, the complicated ε

dependence of the charging free energy in Eq. (3.7) is a direct consequence of the ε

dependence of the solvent charge density in Eq. (3.5).

3.2.2 Gaussian-Smoothed Dipole Model

While Gaussian-smoothed dipole (and multipoles) has been employed to con-

struct force fields in ab initio calculations [61], the motivation of using a Gaussian-

smoothed dipole in the context of solvation is to capture the slowly-varying long-

range component of the potential arsing from an ideal point dipole, allowing us

to focus on the linear component in the static and dynamic responses of dipole

solvation.

A Gaussian-smoothed dipole can be achieved by positioning a Gaussian charge

distribution of magnitude Q centered at the origin, and another Gaussian charge

distribution of magnitude −Q at an infinitesimal distance −dr from the origin.

Recall that a unit Gaussian charge distribution has the same potential as the long-

range component v1(r) in the LMF seperation for the Coulomb potential, Eq. (2.3).

The potentials generated by each individual charge distribution here are Qv1(r) and

−Qv1(r+dr), respectively. Then the potential of the composite Gaussian-smoothed

dipole is given by

vGD (r) = −Q{v1 (r+ dr)− v1 (r)}

= −Qdr · ∇v1 (r) ,

(3.8)

where the product Qdr defines the dipole moment µ, and is kept finite in analogous

40



0 2 4 6 8

r ( )

0.0

0.1

0.2

0.3

vµ
(r

)/
µ
co

sϕ

v µpoint

v
µ
GD, l=2

v
µ
GD, l=3

Figure 3.1: Comparison between the radial part of the point dipole potential vpoint
and the radial part of the Gaussian dipole potential vGD for two different values of
smoothing length l. In the limit as l → 0, vGD reduces to vpoint. With greater l,
vGD becomes more slowly varying.

to the case where a point dipole is constructed from two point charges [62]. Then we

can rewrite the potential of the Gaussian-smoothed dipole with smoothing length l

as

vµGD (r) = −µ · ∇erf (r/l)

r
. (3.9)

Assuming, without loss of generality, that µ = µẑ, Eq. 3.9 then takes the form

vµGD (r, φ) = µ cosφ

[
erf (r/l)

r2
− 2e−r2/l2

√
πlr

]
, (3.10)

which on length scales greater than l, approaches the potential of a point dipole

vpoint(r, φ) =
µ cosφ

r2
. The radial part of vµGD (r, φ) is illustrated in Fig. 3.1. With

the smoothing parameter l chosen on the scale of nearest neighbor correlation of

the system, it is slowly varying with respect to changes in r, thereby avoiding the

41



singularity in the point dipole potential without requiring the shielding, excluded

volume core employed in the Bell model. Moreover, the angular part of vµGD (r, φ) is

transformed to a low degree spherical harmonics Y 0
1 in the Fourier space. Therefore,

we can reasonably consider the potential of a properly constructed Gaussian dipole

to be slowly varying over space and contains only small k components in the k-space,

similar to the behaviors of a Gaussian charge distribution.

Provided that the potential arising from the Gaussian dipole is slowly varying,

it can be considered as an external potential weakly coupled to the aqueous system,

enabling the linear response theory treatment for solvation free energy introduced

in Chapter 2.

Recall in Sec. 2.2, Eq. (2.10) is derived as the linear relation between the

induced solvent charge density ρq(r) and the external potential vQ(r). The Fourier

transform of Eq. (2.10) reads

ρ̂q(k) ≈ −βχ̂qq(k)v̂Q(k). (3.11)

Combining Eq. (3.11) with the k-space version of Poisson’s equation

δv̂(k) =
4π

k2
δρ̂(k), (3.12)

we have

v̂q(k) ≈ −4πβ

k2
χ̂qq(k)v̂Q(k), (3.13)

which relates v̂q(k), the potential due to the induced solvent charge density, to the

external potential. The total electrostatic potential of the system is then obtained

42



by

v̂tot(k) = v̂q(k) + v̂Q(k)

=
(
1− 4πβ

k2
χ̂qq(k)

)
v̂Q(k).

(3.14)

Given that any large k components in v̂Q(k) are attenuated for the slowly-varying

potentials such as those arising Gaussian charge distributions or Gaussian-smoothed

dipoles with sufficiently large l, we consider the long-wavelength limit of k → 0,

which amounts to the macroscopic limit of dielectric screening. Therefore, we have

lim
k→0

v̂tot(k)

v̂Q(k)
= lim

k→0

(
1− 4πβ

k2
χ̂qq(k)

)
=

1

ε
. (3.15)

With Eq. (3.15), we can expand χ̂qq(k) to the second order of k as

χ̂qq(k) ≈ χ̂(0)qq(k) + k2χ̂(2)qq(k), (3.16)

where

χ̂(0)qq(k) = 0 (3.17)

due to the requirement of neutrality, and the second moment of χ̂qq(k) is connected

to the dielectric constant through an generalization of the Stillinger–Lovett moment

conditions [7, 45, 63–65],

4πβχ̂(2)qq(k) =
(
1− 1

ε

)
. (3.18)

The odd moments of χ̂qq(k) vanish due to the nature of the Fourier transform of

χqq(r) in spherical coordinates:

χ̂qq(k) = 4π

∫ ∞

0

sin kr

kr
χqq(r)r2dr. (3.19)

43



By combining Eq. (3.11) and (3.12), we express the the induced solvent charge

density in terms of the solute charge density

ρ̂q(k) ≈ −4πβ

k2
χ̂qq(k)ρ̂Q(k), (3.20)

which, with χ̂qq(k) having been evaluated in Eq. (3.16) – (3.18) for the cases of

perturbation potentials with only small k components, leads to

ρq(r) ≈ −
(
1− 1

ε

)
ρQ(r). (3.21)

To determine the solvation free energy of the Gaussian dipole, we consider the

charging process of a Gaussian dipole in the dielectric continuum by coupling its

potential to a parameter λ varying from 0 to 1, such that the solute charge density

is ρλµGD (r) = λρµGD (r). Eq. (3.21), the linear response theory approximation to the

solvent charge density, is then applied to the Gaussian dipole solvation, yielding

ρqλ (r) ≈ −
(
1− 1

ε

)
ρλµGD (r) , (3.22)

where the charge density distribution of the Gaussian dipole can be obtained by

applying Poisson’s equation to Eq. (3.10),

ρµGD (r) = µ cosφ
2re−r2/l2

π3/2l5
, (3.23)

as illustrated in Fig 3.2.

Presuming that there is no pre-existing solvent charge density in this model,

the charging free energy can be computed with a form analogous to Eq. (3.2),

∆G (µ) =

∫
drvµGD (r)

∫ 1

0

dλρqλ (r). (3.24)

44



4 2 0 2 4

x (Å)
4

2

0

2

4

z (
Å)

l

0.0032

0.0016

0.0000

0.0016

0.0032

Figure 3.2: Charge density of the Gaussian dipole (Eq. (3.23)) scaled by the mag-
nitude of its dipole, ρµGD(r)/µ, in the xz-plane for fixed y = l, with l = 2 Å. White
dashed lines indicate z = ±

√
πl/2, highlighting the length scale of the effective dis-

tance between the positive and negative charge densities comprising the dipole [17].

Inserting Eq. 3.9 and Eq. 3.22 into Eq. 3.24 leads to

∆G (µ) = −
(
1− 1

ε

)
µ2

3
√
2πl3

. (3.25)

In Fig. 3.3, we compare the analytical description of free energy based on Eq. (3.25)

with the free energy computed from MD simulations using the SPC/E water model.

The magnitude of the dipole moment considered here, 8 D, is roughly the change

in the dipole moment upon electronic excitation of coumarin 153, one of the most

widely used probe molecules in experimental solvation dynamics studies [29]. It is

shown that as we increase the smoothing length l (to progressively achieve a slowly-

varying potential of the Gaussian dipole), the simulation results exhibit improved

agreement with the above analytical description. When l reaches 3.0 Å, roughly close

to the nearest neighbor correlation length in the system, the free energy computed

from the MD simulation deviates from the analytical prediction by 8.7 %, suggesting

45



Figure 3.3: Charging free energy of a Gaussian dipole with µ = 8 D and various
values of smoothing length l. Data points correspond to simulation results and
error bars are smaller than the symbol size. The blue solid line is the prediction of
Eq. (3.25). (inset) Comparison between the simulation results and Eq. (3.25) for
l ≥ 3.0 Å on the finer scale.

that the behavior of the system follows linear response to a good approximation.

Using a more conservative smoothing length l = 4.5 Å, the simulation result demon-

strates closer alignment with the analytical prediction drawn from linear response

theory, showing a deviation of only 2.9 %.

The solvation free energy of a Gaussian dipole is an important component of

the total solvation free energy of a molecularly-detailed dipolar molecule. This can

be seen through a judicious rewriting of the solvation free energy [45], wherein a

Gaussian dipole is first inserted into the solvent, then the molecular core is inserted,

and finally the width of the Gaussian is decreased to the point dipole limit or

46



morphed into a more detailed molecular charge distribution. The free energy of this

first step of inserting the Gaussian dipole is expected to be a significant fraction of

the electrostatic contribution to the solvation free energy, and Eq. (3.25) provides an

accurate analytic expression for this quantity. This is especially true for small solute

molecules, where the free energy of inserting an excluded volume core is less than or

comparable to the free energies in Fig. 3.3 in the linear response regime. For large

dipolar solutes, such as charge-neutral proteins, the Gaussian dipole charging free

energy will be an important contribution to the total solvation free energy, although

the relative partitioning of the free energy components will depend sensitively on l,

with l ∼ 1 nm for a small protein [17,66].

3.2.3 Model-Based Dielectric Propensities

Though the Gaussian dipole model and the Bell model both predict the charg-

ing free energy of a dipole to the correct order of magnitude, their descriptions for

free energies exhibit subtle yet profound differences in the manner of relating to the

dielectric permittivity. While the free energy predicted by the Gaussian dipole model

depends on the dielectric permittivity through the factor −
(
1− 1

ε

)
in Eq. (3.25),

the free energy predicted by the Bell model depends on the dielectric permittivity

through − ε−1
2ε+1

in Eq. (3.7).

To further compare the charging free energy of the Gaussian dipole with that

of the Bell model Eq. (3.7), we write Eq. (3.25) in an alternative form

∆G (µ) = −
(
ε− 1

2ε

)
µ2[

(9π/2)1/6 l
]3 . (3.26)

47



The solvation free energies of the two models are first different by a factor of 2ε+1
2ε

.

This factor approaches 3
2
in the limit of vacuum, while it approaches unity for sol-

vents with large dielectric constants such as water. The origin of 2ε+1
2ε

is the different

dielectric continuum treatments applied to the two models. In the Gaussian dipole

model, the linear response theory approximation that we adopted is consistent with

dielectric continuum screening arguments, which leads to the same ε dependence as

models of ion solvation. This is also expected according to the superposition prin-

ciple of electrostatics, provided that the Gaussian dipole is made of two Gaussian

charges. In contrast, as mentioned in Sec. 3.2.1, the ε-dependence of the Bell model

arises from Eq. (3.5), the singular charge distribution of the solvent. Furthermore,

with different orders of multipole moments, Born- and Bell-like models will generally

exhibit different ε-dependencies of the free energy due to the boundary condition

imposed by the excluded volume. In contrast, the Gaussian models do not exhibit

such a multipole-dependent ε-dependence.

Since for each model, the electrostatic interaction energy, monitored in the

processes of solvation dynamics, follows the same ε-dependence of the charging free

energy, the different model-based dielectric propensities in equilibrium also manifest

in the dynamic responses, which will be discussed in detail in the next section.

48



3.3 Linear Response Description for Relaxation Time Scales in Sol-

vation Dynamics

3.3.1 linear Irreversible Processes

Solvation dynamics explores how the dielectric response of the solvent evolves

over time with respect to variations in the solute charge distribution. This could en-

compass scenarios such as introducing an ion or inducing changes in the solute dipole

moment through photoexcitation [29, 30, 67, 68]. The solvation dynamics descrip-

tions developed upon the classic Born and Bell models predict different relaxation

rates of the electrostatic interaction energy for ion and dipole, respectively [57–59].

Here, utilizing the approximation drawn from the linear response theory, we incor-

porate the Gaussian-smoothed charge and dipole models into the phenomenological

nonequilibrium theory. This approach enables us to examine the relation between

the detailed time scales in solvation dynamics and the multipolar property of the

solute charge distributions.

Consider turning on the solute charge distribution instantaneously in an equi-

librated bulk solvent. This introduction of the solute charge distribution, ρQG(r, t),

defines t = 0, and we assume that ρQG(r, t) generates a electrostatic potential that is

sufficiently slowly varying over space, which can be realized by a Gaussian charge

or a Gaussian dipole distribution with sufficiently large smoothing length. We then

focus on how the solvent dynamically responds to the the instantaneous introduc-

tion of the solute charge and evolves to a new equilibrium state. The dynamical

49



response of the solvent is manifest in its charge density ρq (r, t). We can estimate

the time-dependent induced solvent charge density by

ρq (r, t) = −
(
1− 1

εh

)
ρQG (r, t) +

∫ t

−∞
dt′Φ (t− t′) ρQG (r, t′) . (3.27)

The first term in Eq. (3.27) represents an effectively instantaneous response of the

solvent to the introduction of the solute charge, by making analogous to Eq. (3.22)

with ε replaced by the εh. Physically, this instantaneous response results from a

coarse-graining of high frequency modes of the solvent polarization response, such

that the parameter εh contains the effects of these high frequency modes. The con-

nection of εh to the dielectric dispersion ε(ω) and its precise value will be discussed

in more detail below. The second term describes non-Markovian effects due to the

time lag of the polarization of the solvent in response to the time-dependent solute

charge distribution, where Φ (t− t′) is the response function of the considered sys-

tem, which will be determined below. Because the solute distribution will be turned

on instantaneously at t = 0, the integrand is non-zero only for t′ > 0, but we use

this form for generality. Eq. (3.27) follows directly from the framework of the linear

irreversible process [69].

By changing the integration variable from t′ to s = t− t′, we obtain

ρq (r, t) = −
(
1− 1

εh

)
ρQG (r, t) +

∫ ∞

0

ds Φ (s) ρQG (r, t− s). (3.28)

Then the Laplace (one sided Fourier) transform of Eq. (3.28) takes the form

ρq(r, ω) = −
(
1− 1

εh

)
ρQG(r, ω) + Φ(ω)ρQG(r, ω), (3.29)

where Φ(ω) =
∫∞
0

dt exp(−iωt)Φ(t). Equation (3.29) describes the response of the

50



solvent to a charge distribution harmonically oscillating with frequency ω.

Alternatively, this response can be analyzed by generalizing the linear response

approximation Eq. (3.21) to the dynamical domain as

ρq (r, ω) ≈ −
(
1− 1

ε (ω)

)
ρQG (r, ω) , (3.30)

where ε (ω) is the dielectric dispersion, i.e. the frequency-dependent dielectric per-

mittivity. This approximation is based on the condition that ρQG (r, ω) generates a

slowly-varying potential, which holds for Gaussian charge and Gaussian dipole dis-

tributions with sufficiently large smoothing length l. Indeed, recall in Sec. 2.4 we

demonstrated that the solvation dynamics for such a Gaussian charge distribution

accurately follows linear response.

Comparing Eq. (3.29) and Eq. (3.30) yields an expression for the response

function

Φ(ω) =
1

ε (ω)
− 1

εh
. (3.31)

This expression for the response function can be used to described the time depen-

dence of the induced solvent charge density according to

ρq(r, t) = −ρQG(r, t) +

∫ ∞

0

dsρQG(r, t− s)L−1

{
1

ε(ω)

}
, (3.32)

where L−1{f(ω)} denotes the inverse Laplace transform of f(ω). Equation (3.32)

is the linear response theory approximation for ρq(r, t) to a time-dependent solute

charge density, which we expect to be valid for a spatially slowly-varying Gaussian

charge or dipole distribution.

Then we introduce the explicit time dependence of the solute charge density

51



by turning on the solute charge instantaneously at t = 0

ρQG(r, t) = ρQG(r)Θ(t), (3.33)

where Θ(t) is the Heaviside step function, Θ(t) = 0 for t ≤ 0 and Θ(t) = 1 for t > 0.

Under this form of the solute charge density, Eq. (3.32) becomes

ρq(r, t) = −ρQG(r)

[
Θ(t)−

∫ ∞

0

dsΘ(t− s)L−1

{
1

ε(ω)

}]
. (3.34)

In order to further characterize the collective dynamical response of a polar

solvent to the solute charge distribution, we now need to assume an analytic form

for ε(ω). In general, ε(ω) exhibits many timescales encompassing the various types

of motion in a polar solvent. Low frequency modes typically govern the long-time,

collective polarization response that is responsible for dielectric screening of the

solute charge. Here, we are concerned mainly with this low frequency dynamical

response and the longer-time solvent relaxation. High frequency modes include

electronic response (absent in point charge models), inertial and librational motions,

and molecular vibrations, as well as any faster collective polarization fluctuations [28,

29, 70]. In this work, we coarse-grain out these high frequency modes and use a

description of ε(ω) involving a single-Debye relaxation process [71]

ε(ω) = εh +
ε− εh

1− iωτD
, (3.35)

where τD is the Debye relaxation time, and ε is the static dielectric constant of the

solvent. This form of ε(ω) arises from introducing a frequency cutoff, ωh, such that

processes above this frequency are assumed to occur effectively instantaneously. At

the frequency ωh, the Debye process, represented by the second term in Eq. (3.35),

52



is no longer the dominant relaxation channel, and ε(ω) ≈ εh for ω > ωh. Higher

frequency motions, such as additional Debye-like relaxation processes and short-time

inertial Gaussian decays [28, 29, 70, 72–76], are well understood and can be readily

incorporated into more complex forms of ε(ω), but this is beyond the scope of this

thesis.

The inverse Laplace transform of ε−1(ω) can now be evaluated to yield

L−1

{
1

ε(ω)

}
=

εh − ε

τDε2h
e−t/τL +

δ(t)

εh
, (3.36)

where τL is the solvent longitudinal relaxation time defined by

τL =
(εh
ε

)
τD. (3.37)

Finally, inserting Eq. (3.36) into Eq. (3.34), we arrive at an expression for the time-

dependent induced solvent charge density,

ρq(r, t) = −
(
1− 1

ε

)
ρQG(r)Θ(t)−

(
1

ε
− 1

εh

)
ρQG(r)e

−t/τLΘ(t). (3.38)

The first term in Eq. (3.38) corresponds to the equilibrium limit of the solvent

response; the solute charge density is partially screened by the solvent with static

dielectric constant ε. The second term in Eq. (3.38) interpolates between an effective

screening due to the coarse-grained fast-time response at t = 0+ and the equilibrium

limit as t → ∞ on a timescale of τL.

Solvation dynamics is directly described by the nonequilibrium ensemble av-

erage of the time-dependent energy gap, ∆E(t), which is the averaged difference

in the Hamiltonians of the excited and ground states at time t. For the process

considered here, there is no solute prior to insertion, and the energy gap is given by

53



the electrostatic interaction energy,

∆E(t) =

∫
dr

∫
dr′

ρQG(r
′, t)ρq(r, t)

|r− r′|
. (3.39)

Using the solvent charge density given by Eq. (3.38), we obtain

∆E(t) = Θ(t)∆E(Q; l, ε) + e−t/τLΘ(t)
[
∆E(Q; l, εh)−∆E(Q; l, ε)

]
, (3.40)

where the interaction energy, ∆E(Q; l, ε) = 2∆G(Q; l, ε), is equal to twice the charg-

ing free energy defined in the previous sections. As with the charge density, the first

term in Eq. (3.40) corresponds to the equilibrium limit of the interaction energy.

The non-trivial time evolution of the energy is given by the second term, which

exponentially interpolates between the energy change due to fast-time response, via

εh, at t = 0+ and the static, equilibrium interaction energy on the scale of τL.

The time dependence of solvation dynamics is contained in the normalized

nonequilibrium response function, S(t), defined by Eq. (2.1). Within our theory, it

can be obtained by inserting Eq. (3.40) into Eq. (2.1) to yield

S (t) = 1−
1− 1

εh

1− 1
ε

Θ(t) +
1
εh

− 1
ε

1− 1
ε

(e
− t

τL − 1), t ≥ 0. (3.41)

This form of S(t) is equal to unity at t = 0 and decays to zero at long times, as

expected. The second term in Eq. (3.41) corresponds to an effectively instantaneous

decay due to high frequency (ω > ωh) processes. The Debye relaxation is contained

in the third term of Eq. (3.41), which is a single exponential decay over a timescale

of τL. Provided that the response accurately follows linear response, this form of

S(t) is independent of the multipolar property of the charge distribution, e.g. S(t)

is predicted to be the same for a Gaussian charge and a Gaussian dipole.

54



3.3.2 Probing the Dielectric Dispersion of Water with Gaussian-Smoothed

Solute Models

We performed MD simulations of solvation dynamics in the SPC/E model

of water by turning on either a Gaussian charge or a Gaussian dipole at t = 0

and computed S(t) as a nonequilibrium average over many simulation trajectories

following Eq. (2.1). The resulting normalized nonequilibrium response function are

shown in Fig. 3.4 (curves labeled MD).

First, we note that the two simulated S(t) curves exhibit a substantial degree

of overlap, showing that the response to the Gaussian charge and the Gaussian

dipole are essentially the same. This is in agreement with the expectations from the

linear response theory described in Sec. 3.3.1, i.e. the time evolution for solvation

dynamics do not distinguish between a Gaussian charge and a Gaussian dipole,

provided they generate sufficiently slowly-varying potentials.

At very short times, S(t) rapidly decays from one to about 0.25. This fast

decay is known to follow a Gaussian in time and arises from inertial motion in

response to instantaneous turning on of the solute field [29, 70]. Following this

initial response, S(t) exhibits damped oscillations followed by an exponential decay.

The coherent oscillations arise from underdamped, collective solvent motions on

intermediate timescales, wherein solvent molecules initially overpolarize due to the

rapid switching on of the solvent and then overcorrect, with this continuing until the

oscillations are damped on a scale of about 0.5 ps [29]. This damped oscillatory decay

as well as the initial Gaussian decay process are high frequency motions (ω > ωh)

55



Figure 3.4: Comparison between S(t) computed from MD simulations for a Gaussian
charge (black solid line) with Q = 3e0 and l = 5.6 Å, and a Gaussian dipole (red
dashed line) with µ = 8 D and l = 3.0 Å, and S(t) described by Eq. (3.41) (violet
dotted line) with parameters specified in the text. The exponential curve whose
magnitude and rate are determined from the analytical theory can be considered as
the asymptote of the long-time decay in the simulation result. Both plots contain
the same data, the right panel is plotted on a semi-log scale.

and are therefore neglected in the theory for S(t) in Eq. (3.41).

To compare the prediction of Eq. (3.41) with the simulation results, one needs

to plot the analytical prediction with the fitted parameters which physically orig-

inate from the dielectric dispersion in the Debye form, Eq. (3.35), as shown in

Fig. 3.4. The instantaneous drop in S(t) at t = 0+ represented by the second term

in Eq. (3.41) approximates the effects of all responses with frequencies higher than

ωh. This instantaneous drop is then followed by an exponential decay, described by

the third term in Eq. (3.41). The exponential decay captures the lower envelope of

the response at times less than 0.5 ps, indicating that the magnitude of the initial

drop in S(t) in adequately described by εh. At longer times, Eq. (3.41) accurately

captures the long-time decay to equilibrium, with the analytic prediction overlap-

ping with the simulated results, for both the Gaussian charge and dipole. This

56



is emphasized by the examination of S(t) on a logarithmic scale, highlighting the

accuracy of the long-time solvent response predicted by Eq. (3.41).

The coarse-grained initial change described by Eq. (3.41) depends on the rel-

ative magnitudes of εh and ε, and the exponential decay in S(t) is exactly deter-

mined by the complete set of physical parameters governing the dielectric dispersion.

Therefore, fitting S(t) to Eq. (3.41) encodes the information about the dielectric dis-

persion of the solvent. It is indicated that using the dielectric dispersion in the Debye

form with ε ≈ 72, εh ≈ 4, and τD ≈ 8, Eq. (3.41) can capture the response of the

system on different time scales, particularly providing an accurate description of

the long-time decay. The dielectric dispersion obtained from the above fitting pro-

cess agrees well with the results determined for the SPC/E water model in previous

work [77].

We emphasize that our theory for S(t), Eq. (3.41), generally applies to any

time-dependent charge distribution that generates a spatially slowly-varying poten-

tial, like the Gaussian distributions described here. In these models, solute potentials

vary slowly over molecular length scales and induce local and linear responses in the

solvent, such that Eq. (3.30) is valid. For example, Eq. (3.41) correctly predicts

the same time constant τL for the solvation dynamics of both ionic and dipolar

Gaussian-smoothed models, as demonstrated in Fig. 3.4. This longest time scale for

the decay of S(t) is the same as that derived with the classic Born model for ionic

solvation. However, the Bell model leads to a different time constant for solvation

dynamics in response to a dipolar solute composed of a hard core surrounding a

57



Figure 3.5: Analytical predictions derived from our linear response theory-based
model and the traditional models with hard cores. The same parameters for the
single Debye process are used in all the models. Both plots contain the same data,
the right panel is plotted on a semi-log scale.

point dipole,

τBell
L =

(
2εh + 1

2ε+ 1

)
τD. (3.42)

The difference between τL and τBell
L originates from the different ε dependencies

of the free energy discussed in Sec. 3.2.1. In general, the electrostatic boundary

presented by the solute excluded volume leads to multipole-dependent prefactors in

Eq. (3.42). Moreover, the relative magnitudes of the short- and long-time responses

predicted by hard-core models also differ depending on the order of the multipole,

as shown in Fig. 3.5. Therefore, hard-core models lead to non-negligibly different

behaviors for the solvation of different multipoles [57–59], which can be attributed

to the improper application of macroscopic electrostatic approaches to an atomi-

cally sharp boundary. As a result, the hard-core models fail to reflect the intrinsic

dielectric response of the solvent.

58



3.4 Conclusion

Through the examination of the static and dynamic solvation of Gaussian

charges and dipoles in a widely-used model of liquid water, we have demonstrated

that the intrinsic response of water to slowly-varying electrostatic perturbations

is well-described by linear response theory. The use of Gaussian-smoothed charge

distributions removes the singularity associated with point charge models at short

distances, eliminating the need for an excluded volume solute core. The presence of

such cores has obscured the interpretation of solvent response because the harshly

repulsive cores inherently induce non-linear and asymmetric thermodynamic and

structural response in polar solvents [45]. By avoiding these nonlinearities, we are

able to probe and develop theories for the intrinsic dielectric response of water in

the linear regime.

In classic dielectric continuum treatments, the introduction of a solute core

creates a boundary condition that leads to relaxation times containing complicated,

multipole-dependent functions of the dielectric dispersion, which can be consid-

ered as an artifact. In contrast, the solvent response to boundary-less Gaussian

multipolar distributions take on the same dielectric screening form. This leads to

multipole-independent timescales for solvation dynamics, which encode the infor-

mation about the intrinsic solvent response to electrostatic solute perturbations in

the linear response regime, and therefore allow us to probe the dielectric dispersion

of the solvent.

59



Chapter 4: Urea Nucleation in Water I: Exploring Free Energy Land-

scape with Local Molecular Field Theory andWell-tempered

Metadynamics1

4.1 Introduction

In contrast to the processes of solvation dynamics which typically occur on

the picosecond timescale, crystal nucleation takes place on a much greater timescale

and exhibits the rare event nature. When it comes to nucleation from aqueous

solutions, the intermolecular Coulomb interactions can play a significant role in

determining the behavior of the system as well, thus providing a testing ground for

the applicability of LMF theory in complex molecular systems.

Nucleation is the first stage of crystallization, a first-order phase transition in

which structures with long-range correlations form within melts or solutions. It is

of particular importance in fundamental physical and biological scie