J Mol Model (2014) 20:2359
DOI 10.1007/s00894-014-2359-5
ORIGINAL PAPER
Establishing conditions for simulating hydrophobic solutes
in electric fields by molecular dynamics
Effects of the long-range van der Waals treatment on the apparent particle mobility
Zoran Miliˇcevi´c · Siewert J. Marrink ·
Ana-Sunˇcana Smith · David M. Smith
Received: 18 January 2014 / Accepted: 15 June 2014
© Springer-Verlag Berlin Heidelberg 2014
Abstract Despite considerable effort over the last decade,
the interactions between solutes and solvents in the presence
of electric fields have not yet been fully understood. A very
useful manner in which to study these systems is through
the application of molecular dynamics (MD) simulations.
However, a number of MD studies have shown a tremendous
sensitivity of the migration rate of a hydrophobic solute to
the treatment of the long range part of the van der Waals
interactions. While the origin of this sensitivity was never
explained, the mobility is currently regarded as an artifact
of an improper simulation setup. We explain the spread
in observed mobilites by performing extensive molecular
This paper belongs to a Topical Collection on the occasion of Prof.
Tim Clark’s 65th birthday
Electronic supplementary material The online version of this
article (doi:10.1007/s00894-014-2359-5) contains supplementary
material, which is available to authorized users.
Z. Miliˇcevi´c · A.-S. Smith ( ) · D. M. Smith ( )
Cluster of Excellence: Engineering of Advanced Materials,
Friedrich-Alexander University Erlangen-Nuremberg,
Nägelsbachstraße 49b, 91052 Erlangen, Germany
e-mail: smith@physik.fau.de
e-mail: dsmith@irb.hr
Z. Miliˇcevi´c · A.-S. Smith
Institute for Theoretical Physics, Friedrich-Alexander University
Erlangen-Nuremberg, Staudtstraße 7, 91058 Erlangen, Germany
A.-S. Smith · D. M. Smith
Ru ¯der Boškovi´c Institute, Bijeniˇcka cesta 54, 10000 Zagreb,
Croatia
S. J. Marrink
Groningen Biomolecular Sciences and Biotechnology Institute
and Zernike Institute for Advanced Materials, University of
Groningen, Nijenborgh 7, 9747AG Groningen, The Netherlands
dynamics simulations using the GROMACS software package
on a system consisting of a model hydrophobic object
(Lennard-Jones particle) immersed in water both in the
presence and absence of a static electric field. We retrieve
a unidirectional field-induced mobility of the hydrophobic
object when the forces are simply truncated. Careful
analysis of the data shows that, only in the specific case of
truncated forces, a non-zero van der Waals force acts, on
average, on the Lennard-Jones particle. Using the Stokes
law we demonstrate that this force yields quantitative agreement
with the field-induced mobility found within this
setup. In contrast, when the treatment of forces is continuous,
no net force is observed. In this manner, we provide a
simple explanation for the previously controversial reports.
Keywords Molecular dynamics · Electrophoretic
mobility · Van der Waals interactions · GROMACS
Introduction
The hydrophobicity of an interface, droplet or a particle
can be modulated by an external electric field in a process
called electrowetting [1]. Such a change of surface characteristics
allows for the regulation of macroscopic properties
such as adhesion or friction in micro and nano-fluidics
by the electric field [2]. Microfluidic movement [3] based
on electrowetting is being used in an increasing number
of applications, for example, for reflective displays [4, 5].
While the effect of the electric field is reasonably well
understood on the macroscopic scale [1], there are only a
few studies of the origin of electrowetting on the nanometric
scale [6–8]. With the advent of nano-electronics and the
development of micro and nano-devices that can be powered
by electric fields [9], this becomes an increasingly pressing
problem.
/Published online: August 20148
2359, Page 2 of 11 J Mol Model (2014) 20:2359
In the presence of a static electric field, it has been found
that the field exerts torques on non-spherical nano-particles
[10]. This occurs through a coupling between the interfacial
hydrogen bonds and the alignment of the water molecules,
which has also been studied in the context of nanodrops
on graphene surfaces [11] and for water in nanopores
[6, 12]. It is natural to expect that this coupling will also
affect the hydration of spherical nano-particles, and thus
have implications on the self-assembly process, the propulsive
efficiency of nano-devices, and the electrophoretic
mobility [13] of particles.
Recently, an interesting coupling between the electric
field and the solvent organization was suggested, supposedly
leading to net water flux through carbon nanotubes
in the presence of static fields. The latter were directly
imposed [14] or mimicked by an asymmetric distribution of
point charges along the tube [15], while the van der Waals
(vdW) forces were truncated after a certain cutoff distance.
However, in both cases, the flux was found to vanish when
the van der Waals interactions were gradually turned off
[16, 17]. Additionally, the flux was found to be somewhat
sensitive to the neighbor list update frequency, the
commonly used Berendsen thermostat, as well as to the
improper use of the charge groups [17, 18]. However, since
the strongest effect was associated with the treatment of the
long range vdW interactions, the mobility observed with the
truncated forces was argued to be the consequence of the
imprudent implementation of the Lennard-Jones cut scheme
within GROMACS [16, 18, 19].
Similar arguments were evoked after the electrophoretic
mobility of heptane droplets in ion-free water in the presence
of static external fields was reported as a result of an
MD study, again using GROMACS [20]. The mobility was
initially related to the interfacial water structure, involving
static properties such as the water dipolar ordering and the
density profile as well as dynamic properties such as the
viscosity and the slip length [21]. The mechanism was later
disputed by Bonthuis et al [19], who showed that in an
electroneutral dipolar fluid the static electric fields do not
give rise to interfacial flow, even in the presence of dipolar
ordering at the surface. Their conclusion was supported
by the observation that the mobility indeed vanishes when
the long-ranged van der Waals interactions are treated by
the shift approach instead of the default cut(off) approach.
Even though the correctness of the GROMACS package was
questioned several times, quantitative explanation of such
tremendous sensitivity of the nano-fluidic mobility on the
treatment of the van der Waals interactions and on other
simulation parameters has not yet been presented.
In this work we attempt to explain this sensitivity by
constructing a minimal model that could show mobility.
More specifically, as shown in Fig. 1, we perform extensive
MD simulations in GROMACS of a perfectly spherical
Fig. 1 A snapshot from the simulation of a Lennard-Jones particle
(large green sphere) immersed in water
hydrophobic particle in water, both in the presence and
absence of an external electric field.
Methods
Computational details
The standard simulation setup involves a single uncharged
Lennard-Jones (LJ) sphere, with a mass of 50 atomic mass
units and the van der Waals interaction parameters σLJ−O
and εLJ−O of 1.5 nm and 0.8063 kJ mol−1, respectively.
These were deliberately set to create a particle that is
significantly larger than the water molecule (Fig. 1). The
particle was dissolved in a cubic box of edge length 8.9
nm containing 23419 SPC/E [22] water molecules. The
energy of the system was first minimized with the steepest
descent method. Subsequently, a short (0.5-1 ns, NPT) run
was performed followed by an additional equilibration run
(0.5-1 ns, NPT) at nonzero field for the simulations in the
presence of an external electric field. The latter was chosen
to have a strength of 0.6 V nm−1, and was imposed
using tinfoil boundary conditions [23]. During the equilibration,
the pressure was weakly coupled to a target value of
p = 1 bar using the Berendsen barostat [24]. In all cases, the
volume of the system fully relaxed and began to fluctuate
around an equilibrium value.
The production runs (each with a total length of 100 ns)
were conducted in both the NVT and the NPT ensembles,
in the presence and absence of a static electric field
(imposed in the +x direction). For the simulations in the
NVT ensemble, the initial 0.5-1 ns was omitted from the
analysis, to allow the relaxation of the system after fixing its
volume, i.e. density. For all combinations of the ensemble
J Mol Model (2014) 20:2359 Page 3 of 11, 2359
and field strength, the simulations were performed for the
three common treatments of the van der Waals interactions,
namely the cut, switch, and shift treatments.
The standard setup involved simulations with a time step
of 2 fs, the LINCS algorithm [25], particle decomposition,
and a non-bonded list update frequency (NLUF) of 10
steps. Long range electrostatic interactions were accounted
for by the Particle Mesh Ewald (PME) [26] technique. A
reference temperature of 300 K was achieved by employing
the Berendsen algorithm with a time relaxation constant
of 0.1 ps [24]. To ensure good statistics and adequate
spatial and time resolution, the coordinates were saved, with
single precision, every 100 steps. To evaluate the results,
single precision simulations (forces, velocities, and positions)
were compared to runs subject to double precision
evaluations (through all levels of the simulation).
Unless indicated otherwise, simulations were performed
in the GROMACS 4.0.5 simulation package [27].
Additional simulations were executed in GROMACS 3.3.3,
4.5.5, and 4.6.4.
In the study of the effect of the thermostat (NVT ensemble,
NLUF=1, and NLUF=10), the standard Berendsen
weak coupling scheme was replaced with the Nosé-Hoover
[28, 29] thermostat (with a period of oscillations of kinetic
energy set to 0.3 ps), or the advanced velocity rescaling
[30] algorithm (relaxation time of 0.1 ps), as implemented
in the 4.0.5 version of GROMACS. In addition, Langevin
thermostat (friction constant 0.5 ps−1) was implemented
together with the Langevin dynamics [31].
Treatments of the van der Waals interactions
Generally, van der Waals interactions are taken into account
through the well known 12-6 Lennard-Jones potential
(dashed line in Fig. 2a)
VLJ = 4εij
σij
rij
12
−
σij
rij
6
, (1)
where εij is the potential well depth and σij is the distance
between atoms i and j at which VLJ = 0.
An important characteristic of the LJ potential (1) is its
slowly varying dispersion term, that makes the potential
long-ranged. To make computation feasible, the most
common way is to neglect all interactions between atoms
at distances larger than the chosen cutoff distance rc. This
is referred to as the cut treatment. While computationally
very effective, this treatment has a discontinuity both in the
potential and in the force at the cutoff distance (black lines
in Fig. 2).
Because we here focus on the nonbonded interactions
of the LJ particle with water, the cutoff distance is chosen
relatively large. More specifically, we explicitly take into
−0.5
0
LJ potential
cut
switch
shift
1.6 2 2.4 2.8 3.2
r (nm)
−2
−1
0
a)
b)
r1
rc
Fig. 2 The potential (a) and the force (b) for the cut (black lines),
switch (purple lines) and shift (green lines) treatments of the van der
Waals interactions between the oxygen atom of a water molecule and
the LJ particle, when σLJ−O = 1.5 nm, εLJ−O = 0.8063 kJ mol−1.
The dot-dashed lines indicate the boundaries of the transition region,
namely r1 = 2.4 nm and the cutoff radius rc = 2.8 nm
account the interaction of the LJ particle with at least the 4
nearest layers of water
rc σLJ-O + 4σH2O , (2)
where σH2O = 0.31 nm. Hence, the cutoff radius is rounded
to 2.8 nm.
In order to replace the truncated potentials by continuous
ones that also have continuous derivatives, other treatments
of the vdW potential have been proposed [32, 34, 35]. In
the context of the GROMACS package, two such treatments
are extensively used. The first is called shift [34] and it
introduces a function which modifies the potential over its
whole range (0 ≤ r < rc). The second is called switch and
modifies the potential over part of the range (r1 ≤ r < rc).
Fundamentally, there is no difference between the switch
and shift treatments since for r1 = 0 the latter reduces to the
former one. The region r1 ≤ r < rc, which shows the largest
deviations of the potential and/or force from the original,
will be referred to as the transition range.
Within the switch treatment (purple lines in Fig. 2),
the LJ interaction potential (1) is multiplied by a switch
function S(r) defined as
S(r) =
⎧
⎨
⎩
1, if r ≤ r1
1 − 10W3 + 15W4 − 6W5, if r1 < r < rc
0, if r ≥ rc
(3)
2359, Page 4 of 11 J Mol Model (2014) 20:2359
where W = (r −r1)/(rc −r1). This multiplication drives the
potential function towards zero at the cutoff distance. The
van der Waals force in the switch method reads
Fsw
(r) = −∇(VLJS) = FLJS(r) − VLJ∇S(r) , (4)
with FLJ = −∇VLJ(r). The first derivative of the switch
function ∇S(r), and hence the contribution to the total
force, is nonzero only in the transition region. Here, the
force Fsw(r) gains an additional and spurious second
minimum (purple line in Fig. 2b), consequences of which
will become evident at a later stage.
In the shift treatment (green lines in Fig. 2) a function is
added to the LJ potential [32, 34]
V sh
(r) = VLJ(r)−V c
LJ+4εij σ12
ij Arep
− σ6
ij Adis
+V sh
d (r).
(5)
Here, Vc
LJ is the value of the LJ potential at the cutoff
distance rc, while Arep and Adis are the repulsion and
dispersion corrections of the LJ potential depending only on
the values of r1 and rc. The term V sh
d is nonvanishing only
for r1 < r < rc and equals
V sh
d (r) = −
4
3
εij σ12
ij K
rep
1 + σ6
ij Kdis
1 (r − r1)3
− εij σ12
ij K
rep
2 + σ6
ij Kdis
2 (r − r1)4
,
(6)
where K
rep
1 , Kdis
1 , K
rep
2 and Kdis
2 are constants depending
only on the values of r1 and rc and the nature of the interaction,
i.e. repulsive or attractive (see Supporting Information
for full expressions). The resulting potential V sh(r) deviates
only slightly from VLJ(r) (green line in Fig. 2a). Importantly,
the resulting force Fsh(r) deviates from FLJ(r) only
in the transition region, where it smoothly approaches zero
at rc due to the contribution Fsh
d (r) = −∇V sh
d (r) (green line
in Fig. 2b).
Results and discussion
Water at the interface of the LJ particle
The static structural properties of water around a model
hydrophobic object, i.e. the Lennard-Jones particle, can be
quantified by the radial distribution function (rdf) between
the particle and the oxygen atom of the water molecule gLJ-O
[33, 36–38, 41].
In Fig. 3 we show the rdfs gLJ-O for all considered treatments
of the vdW interactions in the absence of an external
electric field. Although we present data obtained in the NVT
ensemble, identical results are obtained in the NPT ensemble.
The common main features are three well-defined
maxima and minima, the latter indicating the edges of the
respective hydration shells. Strikingly, a careful inspection
1.6 2 2.4 2.8
r (nm)
0.92
1
1.08
1.16
gLJ-O
(r)
cut switch shift
2.4 2.8
1
1.01
1.56 1.6
1.12
1.14
1.16
1.72 1.76
0.94
0.96
1.84 1.92
1.02
1.03
1.04
1st
maximum
1st
minimum
2nd
maximum
Fig. 3 Comparison of the radial distribution function g(r) between the
LJ particle and water oxygen for the cut, switch and shift treatments
of the van der Waals interactions. Inset the behavior of g(r) at long
distances, r ≥ 2.1 nm. Bottom row the local behavior of g(r) at the
positions of 1st and 2nd maxima, and 1st minimum, respectively. The
results are for the simulations performed in the NVT ensemble in the
absence of an electric field
of the transition region (see inset of Fig. 3) reveals an additional
and unexpected structuring in the case of the switch
treatment compared to the cut and the shift treatment within
the transition region. A clear indication that these correlations,
although small, are unphysical is the appearance of
more negative correlations (g < 1) at the position of the
fourth minimum than at the third minimum. Furthermore,
even though the vdW force is identical for r < r1 in all
three approaches, the perturbation imposed in the transition
region propagates towards the hydrophobic object and
manifests in subtle changes of water structure closer to the
hydrophobic object, as evidenced in the small graphs in the
lower row in Fig. 3 representing gLJ-O(r) around the 1st
and 2nd maxima, as well as the 1st minimum. The same
increase of correlations observed for the switch treatment in
gLJ-O is apparent in the rdf between the LJ particle and the
water hydrogens gLJ-H (some examples are given in the Supporting
Information, Fig. S1a). Furthermore, this increase is
independent of the version of GROMACS and the precision
at which the simulation is executed.
It is important to stress that similar spurious correlations
were not found in any of the water-water rdfs (i.e. gO-O,
gH-H, and gO-H), supposedly because the employed cutoffs
are too large to have impact on the correlations of the much
smaller solvent molecules. However, even with transition
regions at much shorter distances (0.6 − 0.9 nm, 0.7 − 1.0
nm, 0.9 − 1.2 nm, 1.5 − 1.8 nm), the described increase of
J Mol Model (2014) 20:2359 Page 5 of 11, 2359
correlations was not found in the rdfs of the water molecules
in simulations of a pure water box and the switch treatment.
This is presumably because at these distances the dominant
interactions between water molecules are of Coulomb type.
In the presence of the electric field, the analysis of
the water structure around a hydrophobic object has not
yet been performed systematically for different vdW treatments.
The difficulty is the breaking of the isotropy due
to a specifically imposed field orientation. Consequently,
averaging over two spatial coordinates respecting azimuthal
symmetry must be performed to fully resolve the distribution
of water at the interface. Unfortunately, such an
analysis would require sampling statistics that are currently
out of reach. We thus analyze the radially averaged distribution
function which still contains the relevant information
(Fig. S1b). We find that when vdW interactions are treated
with the switch function, the same increase of structuring
as in the absence of the field can be seen. Such an observation
is in agreement with previous results on the treatment
of electrostatic interactions with an atom-based switch
function [42]. Therein, it was shown that such an electrostatic
switch treatment induces artificial dipolar ordering of
water in the transition region [35].
It is clear that different treatments of vdW interactions
indeed induce subtle changes of the water structure around
a hydrophobic object. However, it is not clear if these
small differences can be related to the potentially different
−40
−30
−20
−10
0
10
Displacement(nm)
cut
switch
shift
40 60 80
−40
−30
−20
−10
0
40 60 80
t (ns)
0 20 0 20
E0
NVT ensemble
NPT ensemble
Fig. 4 The displacement of the LJ particle in the absence (left) and the
presence of static electric fields (right). In the absence of electric fields
the average displacement over three Cartesian directions is shown. In
the presence of electric fields the displacement in the field direction
(imposed in the +x direction) is shown. In the upper and lower row
are data for the NVT and the NPT ensemble, respectively. The total
simulation time for each setup is 100 ns. The dashed black lines are
the linear fits of the displacement to illustrate the unidirectionality of
the movement of the Lennard-Jones particle
mobilities of LJ particles in water. To see these effects, we
further study the dynamics of the system.
Observed time-dependent displacement of the LJ particle
We first investigate the time-dependent displacement of
the LJ particle both in the presence and absence of a static
electric field (Fig. 4). At zero field, the average displacement,
= ( x + y + z)/3, is fully comparable to a
trajectory along a particular coordinate. In the presence of
the electric field, we focus on the displacement along the x
axis which coincides with the direction of the field.
In Fig. 4, we show the results for both the NVT and
the NPT ensembles and all three treatments of vdW interactions.
Independently of the ensemble, the majority of
the displacement curves show a strong resemblance to
trajectories expected from a particle exhibiting Brownian
motion. The switch treatment provides an analogous
behavior to the shift, despite the slight differences in the
static organization of interfacial waters. On the time scale
of 100 ns, the displacements observed with these two treatments
in the presence of the field are indistinguishable
from those at zero field. The exception is the superficially
small preference to move in a direction opposite to the
field, which is an effect that we could not prove to be
statistically relevant at this stage. However, a stark contrast
can be seen in the displacements under the field when
using the truncated vdW forces. These curves show a linear
drift superimposed on the Brownian displacement, which
is suggestive of a net force acting on the particle. We find
field-induced velocities of the LJ particle of −0.45 m s−1
(NVT) and −0.41 m s−1 (NPT), obtained from linear fits to
the respective data sets (dashed lines).
Comparable velocities were found in several previous
studies with the cut treatment [15, 17–19, 21] on similar
systems. The latter was declared spurious after no net movement
(or water flux) was found on comparable time scales
when the shift method was employed with the field [17, 19].
Actually, a non-zero flux of water molecules through a carbon
nanotube [18] was observed for the cut treatment even
in the absence of the electric field. Since we find no particular
displacement of the LJ particle when E = 0, the
reported flux seems to be related to the particularity of the
investigated system.
Apart from the cut-off treatment, other parameters such
as the time step, the neighbor list update frequency (NLUF)
or the Berendsen thermostat (BT) were suggested to induce
unphysical motion of water around a hydrophobic object
in GROMACS, although these results were based on
relatively short simulations [17]. We test these parameters
(BT, NLUF) in a set of 100 ns long simulations
(cut treatment, NVT ensemble, and E0 = 0.6 V nm−1), and
show the results in Fig. 5, alongside the reference setup
(BT, NLUF = 10). A constant negative drift velocity is
2359, Page 6 of 11 J Mol Model (2014) 20:2359
40 60 80
t (ns)
−40
−30
−20
−10
0
Displacement(nm)
0 20
Fig. 5 The displacement of the LJ particle in the direction of the
imposed static electric field, Ex = 0.6 V nm−1. The simulations are
performed in the NVT ensemble and the vdW interactions are represented
by the cut scheme. Legend BT - Berendsen thermostat, VR velocity
rescale, NLUF - neighbor list update frequency (# timesteps)
observed for each setup, while the difference between the
curves is within the expected statistical noise. The insensitivity
of the drift motion on the NLUF indicates that the
water molecules positioned around rc, that actually contribute
to the update of the neighbor list of the LJ particle,
do not cause the drift. Thus, the mechanism that drives the
net movement of the LJ particle in a direction opposite to
the field is most likely related to the solute-solvent interface
effects. Additionally, we obtained very similar displacements
with the advanced velocity rescaling (VR) algorithm
[30] (−0.39 m s−1 for NLUF = 10 and −0.46 m s−1 for
NLUF = 1) and the Berendsen algorithm (−0.45 m s−1
for NLUF = 10 and −0.35 m s−1 for NLUF = 1). Hence,
the observed field induced mobility can not be assigned to
the imperfections of the Berendsen thermostat, which does
not rigorously reproduce a correct thermodynamic ensemble
[43].
In addition, we performed simulations of 20 ns, to test
the influence of the system size, the time step, and the field
strength on the observed magnitude of the field-induced
velocity of the LJ particle, in the NVT ensemble using the
cut treatment, the Berendsen thermostat, and NLUF = 10. In
the system with NH2O = 63596, the observed drift velocity
of the LJ particle remained the same as in the smaller
system. With the smaller time step the drift was somewhat
larger, while with the decrease of the field strength to 0.3 V
nm−1 the drift was significantly reduced.
Brownian diffusion of the LJ particle
To elucidate the origin of the observed drift, we first analyze
the underlying Brownian displacement and determine the
diffusion coefficient D of the LJ particle. At zero field, the
diffusion coefficient is obtained from the slope of the three
dimensional mean square displacement as D = (rLJ(t) −
rLJ(0))2 /(6t), where rLJ(t)−rLJ(0) is the distance traveled
by the LJ particle in time t. In the presence of the field, we
calculate the diffusion coefficient in the field direction as
Dx = (xLJ(t) − xLJ(0))2 /(2t), where xLJ(t) is the timedependent
position of the LJ particle along the x axis. In
order to eliminate effects of the drift, and assure linearity,
only the first 50 ps of the mean square displacement curve
were used in all cases.
Table 1 summarizes the obtained diffusion coefficients
of the LJ particle. At zero field, within the statistical uncertainty,
D is found independent of the vdW interaction
treatment and the thermodynamic ensemble employed. The
same applies in the presence of the electric field with the
diffusion coefficient in the field direction possibly slightly
enhanced. Similarly, the diffusion coefficient in the directions
perpendicular to the field (data not shown) remains
unaffected by the treatment of the vdW interactions.
It is worth noting that the measured diffusion coefficients
were not sensitive to the version of GROMACS and most
of the results fall within the statistical significance of one
another. Likewise, the difference between the NVT and
NPT ensembles borders with the statistical accuracy. This is
presumably the consequence of the employed large system
size, when the averages of macroscopic observables become
identical in the two ensembles.
The net vdW force on the LJ particle
The linear displacement with the cut approach is indicative
of a static net force acting on the particle. As the LJ particle
carries no charge, it is difficult to envisage a simple electrostatic
force. The only remaining possibility is a force
originating in vdW interactions, which is further supported
by the sensitivity of the mobility to the vdW treatment.
We thus calculate the instantaneous total van der Waals
force FvdW(t) in each Cartesian direction exerted on the LJ
particle by the solvent, by adding the contributions from all
Table 1 The diffusion coefficients (in 10−6 cm2 s−1) of the LennardJones
particle for different treatments of the vdW interactions in the
presence and absence of electric fields for both the NVT and NPT
ensembles. In the presence of the field, the diffusion coefficient in the
direction of the field is presented
Treatment E0 (V nm−1) NVT NPT
CUT 2.64 ± 0.06 2.58 ± 0.01
SWITCH 0.0 2.61 ± 0.03 2.66 ± 0.05
SHIFT 2.64 ± 0.03 2.58 ± 0.01
CUT 2.70 ± 0.05 2.87 ± 0.05
SWITCH 0.6 2.64 ± 0.01 2.67 ± 0.01
SHIFT 2.62 ± 0.05 2.64 ± 0.03
J Mol Model (2014) 20:2359 Page 7 of 11, 2359
water molecules found within rc. The net vdW force is taken
to be the time average of instantaneous forces
FvdW
i =
1
T
T
j=1
Nw
k=1
FvdW
i (rkj ; rkj < rc) , i ≡ {x, y, z}
(7)
where T is the total number of frames (system configurations)
analyzed and equals 500 000 in each case, Nw is the
total number of water molecules present in the system, and
rkj is the relative position of the oxygen atom (position of
the vdW interaction site) of water molecule k with respect
to the position of the LJ particle at time frame j. The specificity
of the van der Waals treatment (cut, switch, shift)
employed in each particular setup is taken into account by
using an appropriate expression for the force.
Table 2 shows the calculated net vdW forces in each
Cartesian direction on the LJ particle. In the absence of the
external electric field, independent of the treatment of the
vdW interactions and the thermodynamic ensemble, the net
force in each Cartesian direction is effectively zero. This
is consistent with the observed lack of net drift of the LJ
particle. Similarly, with the switch and shift treatment at
E0 = 0.6 V nm−1 the net force remains effectively zero in
all directions, as well as for the cut treatment in the off-field
directions.
Strikingly, in the field direction we find a significant
net force of ≈ −6.6 pN acting on the LJ particle. This
force is quickly converging and exhibits a standard deviation
of about 0.3 pN (see Fig. S2). While this force appears
large, it is of a similar magnitude to that exerted by a motor
molecule on vesicular cargo. In the aqueous context, this
force is about the same as the one acting between the LJ
particle and a single water molecule at a distance of about
1.6 nm. Importantly, the magnitude of this force is not
related to the version of GROMACS, and is only associated
with the cut treatment of the vdW interactions. This is
evidenced by the mean force of about −6.6 pN determined
in simulations with the standard setup in all studied versions
of GROMACS (i.e. 3.3.3, 4.5.5 and 4.6.4). In addition, test
simulations performed with double precision did not affect
the result. The observed force with the cut treatment was
again −6.5±0.1 pN, while the force was vanishing with the
shift and switch treatments.
In a more detailed analysis, we find that this force is
generated by water molecules belonging to the first two
hydration shells while the contribution from waters in the
transition region conforms to a Gaussian distribution with
zero mean. This is true for all setups, including the ones with
the switch treatment.
Consistency between the evaluated vdW force
and the observed drift
If the evaluated force observed with the cut treatment is
consistent with the observed drift velocity, the Stokes law
should apply
FvdW
i = ξi · νi . (8)
Here νi is the drift velocity in the direction i of the net force
acting on the particle. The proportionality constant ξ is the
friction coefficient and is simply related to the diffusion
constant by the Einstein relation ξi = kBT /Di (kB is the
Boltzmann constant).
Using Eq. 8, with Dx from Table 1, and FvdW
x from
Table 2, for the cut treatment under the field, we predict a
drift velocity of −0.43 m s−1 and −0.46 m s−1 for the standard
setup of parameters in the NVT and the NPT ensemble,
respectively. This is in exceptional agreement with the fieldinduced
drift observed in the simulations of −0.45 m s−1
(NVT) and −0.41 m s−1 (NPT). The small difference
between the predicted and the observed drift velocity can be
attributed to the intrinsically present thermal noise resulting
from the Brownian motion of the LJ particle. This agreement
clearly demonstrates that the observed field-induced
mobility of a hydrophobic object is a result of a fluctuating
Table 2 The net van der Waals force on the Lennard-Jones particle in each Cartesian direction for different treatments of the vdW interactions in
the presence and absence of static electric fields imposed in the +x direction
Treatment E0 (V nm−1) NVT ensemble NPT ensemble
FvdW
x (pN) FvdW
y (pN) FvdW
z (pN) FvdW
x (pN) FvdW
y (pN) FvdW
z (pN)
CUT −0.02 ± 0.07 −0.02 ± 0.06 0.00 ± 0.02 −0.01 ± 0.07 −0.02 ± 0.10 0.01 ± 0.05
SWITCH 0.0 −0.05 ± 0.06 −0.01 ± 0.15 0.05 ± 0.08 0.03 ± 0.04 0.02 ± 0.15 0.01 ± 0.13
SHIFT −0.02 ± 0.06 0.00 ± 0.09 0.01 ± 0.11 0.01 ± 0.08 0.07 ± 0.11 0.00 ± 0.08
CUT −6.59 ± 0.04 0.07 ± 0.09 0.02 ± 0.07 −6.62 ± 0.06 −0.03 ± 0.15 0.04 ± 0.11
SWITCH 0.6 −0.03 ± 0.09 −0.04 ± 0.11 0.00 ± 0.13 0.03 ± 0.08 −0.02 ± 0.02 0.05 ± 0.13
SHIFT 0.03 ± 0.06 0.01 ± 0.07 −0.02 ± 0.08 0.03 ± 0.05 −0.01 ± 0.13 0.02 ± 0.10
2359, Page 8 of 11 J Mol Model (2014) 20:2359
but, on average, non-zero force that is produced by truncating
the vdW forces. This force is directly related to the field
strength, as evidenced by its decrease to −4.99 ± 0.03 pN
at the field strength of 0.3 V nm−1.
Consistent with the observed drifts discussed in previous
sections, we find the average net vdW force on the LJ
particle in the field direction for the cut treatment to be independent
of the chosen thermostat and NLUF employed. The
forces obtained for the test set of simulations in the NVT
ensemble (see Fig. 5) are: −6.7 ± 0.1 (VR, NLUF = 10),
−6.61 ± 0.07 (VR, NLUF = 1), and −6.58 ± 0.05 pN (BT,
NLUF = 1). This is in excellent agreement with Fcut
x =
−6.59±0.04 pN obtained for the standard setup (BT, NLUF
= 10). Naturally, the calculated νx arising from the FvdW
x
(between −0.43 m s−1 and −0.46 m s−1) agrees well with
the observed induced field velocities (between −0.35 m s−1
and −0.46 m s−1). With the increase of the system size to
NH2O = 63596 we obtain Fcut
x = −6.72 ± 0.03 pN which
is, expectedly, almost the same as in the smaller system. The
smaller time step of 1 fs slightly increases the net force to
−7.0 ± 0.3 pN, which is in agreement with the somewhat
larger drift velocity observed. Since the diffusion coefficient
Dx is effectively independent of the simulation setup (data
not shown), this data demonstrates that the average net vdW
force found on the hydrophobic object is exclusively related
to the cut treatment of the vdW interactions.
The effects of the thermostat and the “flying ice cube”
effect
In a number of experimental realizations of electrophoretic
experiments, the total number of particles is preserved and
the overall system is canonical. Consequently, only the
average total energy is preserved. One can introduce this to
molecular dynamics by the application of appropriate thermostats,
which modify the Newtonian MD scheme such that
a statistical ensemble is generated at a constant temperature.
However, it is well documented that inappropriate choices
of thermostats may significantly affect the fluctuations in
the system and, in some cases, induce energy drifts caused
by the accumulation of numerical errors [39, 40].
To study these effects in our system, we have performed
a set of simulations in a standard setup with different thermostats
at 300 K. We sample the instantaneous velocity
νx(LJ) along the direction of the imposed field for 10 ns.
Since the system is expected to be canonical, this distribution
should be a Gaussian identical to the one-dimensional
Maxwell-Boltzmann distribution of velocities of a colloid
with a mass associated with the LJ particle, generated at 300
K (symbols in Fig. 6a).
Since it is expected that the canonical nature of the
system is best reproduced by a Langevin thermostat, we perform
Langevin dynamics simulations [31]. In this stochastic
but ergodic scheme, each particle is coupled to a local heat
bath, which removes heat trapped in localized modes. While
the momentum transfer in such simulations is destroyed,
making the diffusion coefficients inaccessible with this thermostat,
the canonical distribution is obtained accurately, as
evidenced in Fig. 6.
In the context of Newtonian dynamics, an attractive and
elegant thermostat is that of Nosé-Hoover [28, 29]. It is
deterministic, time reversible and in principle canonical
(extra degree of freedom acts as thermal reservoir), but
sometimes computationally demanding. However, when
the system is poorly ergodic, this approach may become
difficult and a chain of thermostats may have to be imposed.
However, this does not seem to be the case in the current
simulations, and an excellent agreement with the canonical
−600 −300 0 300 600
vx
(LJ) (m s−1
)
0
0.5
1
1.5
2
Probabilitydensity(10−3
sm−1
)
BT
VR
NH
LA
MB distribution
−600 −300 0 300 600
vx
(LJ) (m s−1
)
a) b)
BT-shift
BT-cut
LA-cut
LA-shift
Fig. 6 (a) Comparison of the distribution of the instantaneous velocity
(saved at every step for 10 ns) of the LJ particle in the direction
of the imposed static electric field, Ex = 0.6 V nm−1 with the shift
treatment. The simulations were performed with the Berendsen thermostat
(BT), advanced velocity rescaling (VR), Nosé-Hoover (NH),
and by Langevin dynamics (LA). The expected Maxwell-Boltzmann
distribution (MB) is shown with symbols. The same analysis for the cut
treatment is presented in Fig. S3. (b) Comparison of the cut (dashed
lines) and shift (full lines) treatments with the Berendsen (black) and
the Langevin thermostats (orange)
J Mol Model (2014) 20:2359 Page 9 of 11, 2359
distribution is obtained (Fig. 6a). An equally good agreement
with the Maxwell-Bolzmann distribution is obtained
with the advanced velocity rescaling method [30], within
which the target temperature is drawn from the canonical
probability distribution, to produce the correct ensemble.
We compare the performance of these canonical thermostats
to the commonly used Berendsen weak coupling
scheme [24]. The Berendsen approach gained popularity
because of its robustness, speed and easy implementation.
The problem with this approach is that it is not timereversible
or deterministic and the resulting distributions
do not strictly correspond to any ensemble. In practice,
however, the deviations from the canonical distribution
are relatively small and decrease with increasing system
size. Specifically in our system, the distribution generated
with this treatment is within the noise of distributions
generated with more formal thermostats for the shift treatment
(Fig. 6a). For the cut treatment the Nosé-Hoover, the
advanced velocity rescaling and the Berendsen thermostats
show small deviations around the maximum from the distribution
obtained with the Langevin thermostat, the latter
coinciding with the distributions obtained with the shift
treatment (Fig. 6b).
The use of the Berendsen thermostat, as well as the
imposition of truncation schemes for the vdW forces is
sometimes associated with localized and unwanted correlated
motions. Spurious drifts may appear due to the transfer
of energy from fast to slow degrees of freedom, which is
known as the “flying ice cube” effect. There are several
steps which help in avoiding this effect, including the use
of the here chosen Verlet algorithm to propagate the system.
Since this algorithm is time reversible, it is associated with
little or no long-term energy drain. Furthermore, it is important
to remove all motion of the center of mass of the system
(translational and rotational) which we perform in each
molecular dynamics step [40]. Beside making relevant
choices in setting the simulations, one can test for the
“flying ice cube” effect a posteriori, by carefully analyzing
the entirety of the data, as discussed below.
Several arguments can be used to support our interpretation
that the drift velocity is associated with the net
force acting on the particle, and not with the “flying ice
cube” effect. The first is that the our system reproduces the
canonical distribution very well. Its size is such that the thermodynamic
limit is approached in the sense that differences
between system averages in the NVT and the NPT ensembles
are basically within the statistical accuracy. Furthermore,
we obtain the same results (net force on the particle as
a function of the treatment of the van der Waals forces) independently
of the thermostat. Specifically, in a set of 10 ns
long simulations, we obtain a net average force in the direction
of applied electric field of −6.7 ± 0.2 pN, −6.5 ± 0.1
pN, and −6.6 ± 0.3 pN with the Berendsen, the advanced
−600 −400 −200 0 200 400 600
vx
(LJ) (m s−1
)
0.5
1
1.5
2
Probabilitydensity(10−3
sm−1
)
0−20 ns
20−40 ns
40−60 ns
60−80 ns
80−100 ns
–
Fig. 7 The distribution of the velocity of the LJ particle along the
direction of the field ¯vx(LJ), calculated from the displacement of the
LJ particle with the resolution of 0.2 ps. The simulation was performed
using the default setup (NVT, BT, NLUF = 10), with the vdW interactions
treated using the shift approach. The distribution was generated
over 5 successive intervals of 20 ns. The comparison of the results
shows no time evolution of the velocity distribution. The same analysis
for the cut treatment is presented in Fig. S4
velocity rescaling and the Nosé-Hoover thermostats, respectively.
Even the Langevin dynamics yields an average force
of −6.9 ± 0.3 pN with the cut treatment. On the other hand
the shift treatment always leads to zero force, independently
of the thermostat and the integration algorithm. Furthermore,
the obtained drift is fully consistent with the observed
force, which only occurs with the cut treatment.
While we are able to confidently associate the spurious
motion of the LJ particle in the electric field to the cut treatment
of the vdW interaction, slow energy transfer could
still be possible on even larger time scales. To exclude this,
we test for the slow draining of the kinetic energy into the
translational degrees of freedom of the particle, by checking
its velocity distribution generated in intervals of 20 ns over
100 ns (original protocol, Fig. 7). Here, the velocity was calculated
from the actual displacement of the LJ particle with
a resolution of 0.2 ps. Importantly, the velocity distribution
is fully independent of time, with any of the treatment of
the van der Waals interactions (see Fig. 7 for the shift treatment
and Fig. S4 for the cut treatment), which we believe,
convincingly excludes the existence of drifts associated with
“flying ice cube” effect in our simulations.
Conclusions
Using a series of molecular dynamics simulations in
GROMACS, we have investigated the behavior of a solvated
nanosized Lennard-Jones particle in the presence of a static
electric field, as a function of the treatment of the longrange
vdW interactions. Even in the absence of the field, we
showed that the different treatments result in subtle changes
of the distribution of water molecules at the interface with
the LJ particle. These changes are not necessarily correct,
2359, Page 10 of 11 J Mol Model (2014) 20:2359
as exemplified by the excessive correlations in the transition
region when the switch treatment is used. The observation
of a fluctuating but non-zero vdW force on the LJ particle
occurring in static electric fields with the cut treatment can
be considered in a similar way. Indeed, it is likely that similar
forces were responsible for the negative electrophoretic
mobilities of hydrophobic objects, and water fluxes through
carbon nanotubes, which were reported in previous studies.
While the accuracy of GROMACS was brought into
question, it seems to be now clear that these effects arise
merely from the truncation of the forces, while the spread
of the previously reported results may also be associated
with the relatively small system sizes and simulation times.
This is further supported by the fact, neither the artifacts
in the static radial distribution function with the switch
treatment, nor the measured forces with the cut treatment
depend on the used version of GROMACS. Analogously,
the same reasoning can be used to rationalize why the shift
implementation for the vdW interactions apparently provides
the most accurate results. Despite these clarifications,
the subtleties in the organization of water in the first two
hydration shells, which give rise to these effects, have yet to
be properly characterized.
Acknowledgments We thank V. Knecht, K. Mecke und U. Felderhof
for stimulating discussions. ASS and DMS acknowledge financial
support from the Cluster of Excellence: Engineering of Advanced
Materials, Erlangen, Germany. ZM acknowledges the partial financial
support from BAYHOST.
References
1. Shamai R, Andelman D, Berge B, Hayes R (2008) Water,
electricity, and between . . . On electrowetting and its applications.
Soft Matter 4:38–45
2. Robinson L, Hentzell A, Robinson ND, Isaksson J, Berggren M
(2006) Electrochemical wettability switches gate aqueous liquids
in microfluidic systems. Lab Chip 6:1277–1278
3. Rotenberg B, Pagonabarraga I (2013) Electrokinetics: insights
from simulation on the microscopic scale. Mol Phys 111:827–842
4. Hayes RA, Feenstra BJ (2003) Video-speed electronic paper based
on electrowetting. Nature 425:383–385
5. Pollack MG, Fair RB, Shenderov AD (2000) Electrowetting-based
actuation of liquid droplets for microfluidic applications. Appl
Phys Lett 77:1725–1726
6. Bratko D, Daub CD, Luzar A (2009) Water-mediated ordering of
nanoparticles in an electric field. Faraday Discuss 141:55–66
7. Dzubiella J, Hansen JP (2005) Electric-field-controlled water
and ion permeation of a hydrophobic nanopore. J Chem Phys
122:234,706–14
8. Vaitheeswaran S, Yin H, Rasaiah JC (2005) Water between plates
in the presence of an electric field in an open system. J Phys Chem
B 109:6629–6635
9. Chang ST, Paunov VN, Petsev DN, Velev OD (2007) Remotely
powered self-propelling particles and micropumps based on
miniature diodes. Nat Mater 6:235–240
10. Daub CD, Bratko D, Ali T, Luzar A (2009) Microscopic dynamics
of the orientation of a hydrated nanoparticle in an electric field.
Phys Rev Lett 103:207,801–4
11. Daub CD, Bratko D, Leung K, Luzar A (2007) Electrowetting at
the nanoscale. J Phys Chem C 111:505–509
12. Bratko D, Daub CD, Leung K, Luzar A (2007) Effect of field
direction on electrowetting in a nanopore. J Am Chem Soc
129:2504–2510
13. O’Brien RW, Beattie JK, Djerdjev AM (2014) The electrophoretic
mobility of an uncharged particle. J Colloid Interface Sci 420:70–
73
14. Joseph S, Aluru NR (2008) Pumping of confined water in carbon
nanotubes by rotation-translation coupling. Phys Rev Lett
101:064,502–4
15. Gong X, Li JY, Lu H, Wan R, Li JC, Hu J, Fang H (2007) A
charge-driven molecular water pump. Nat Nanotechnol 2:709–
712
16. Bonthuis DJ, Falk K, Kaplan CN, Horinek D, Berker AN, Bocquet
L, Netz RR (2010) Comment on Pumping of confined water in
carbon nanotubes by rotation-translation coupling. Phys Rev Lett
105:209,401–1
17. Wong-ekkabut J, Miettinen MS, Dias C, Karttunen M (2010)
Static charges cannot drive a continuous flow of water molecules
through a carbon nanotube. Nat Nanotechnol 5:555–557
18. Bonthuis DJ, Rinne KF, Falk K, Kaplan CN, Horinek D, Berker
AN, Bocquet L, Netz RR (2011) Theory and simulations of water
flow through carbon nanotubes: prospects and pitfalls. J Phys:
Condens Matter 23:184,110–10
19. Bonthuis DJ, Horinek D, Bocquet L, Netz RR (2009) Electrohydraulic
power conversion in planar nanochannels. Phys Rev Lett
103:144,503–4
20. Knecht V, Levine ZA, Vernier PT (2010) Electrophoresis of
neutral oil in water. J Colloid Interface Sci 352:223–231
21. Knecht V, Risselada HJ, Mark AE, Marrink SJ (2008)
Electrophoretic mobility does not always reflect the charge on an
oil droplet. J Colloid Interface Sci 318:477–486
22. Berendsen HJC, Grigera JR, Straatsma TP (1987) The missing
term in effective pair potentials. J Phys Chem 91:6269–6271
23. Hünenberger P (1999) Simulation and theory of electrostatic
interactions in solution: computational chemistry, biophysics, and
aqueous solutions. Am Inst Phys:17–83
24. Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak
JR (1984) Molecular-dynamics with coupling to an external bath.
J Chem Phys 81:3684–3690
25. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM (1997) LINCS:
A linear constraint solver for molecular simulations. J Comput
Chem 18:1463–1472
26. Darden T, York D, Pedersen L (1993) Particle mesh Ewald: An
N · log(N) method for Ewald sums in large systems. J Chem Phys
98:10089–10092
27. van der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE,
Berendsen HJC (2005) Gromacs: Fast, flexible and free. J Comp
Chem 26:1701–1718
28. Hoover WG (1985) Canonical dynamics: equilibrium phase-space
distributions. Phys Rev A 31:1695–1697
29. Nosé S (1984) A molecular dynamics method for simulations in
the canonical ensemble. Mol Phys 52:255–268
30. Bussi G, Donadio D, Parrinello M (2007) Canonical sampling
through velocity rescaling. J Chem Phys 126:014,101–7
31. van Gunsteren WF, Berendsen HJC (1988) A leap-frog algorithm
for stochastic dynamics. Mol Sim 1:173–185
32. Baron R, Trzesniak D, de Vries AH, Elsener A, Marrink SJ, van
Gunsteren WF (2007) Comparison of thermodynamic properties
of coarse-grained and atomic-level simulation models. Chem Phys
Chem 8:452–461
J Mol Model (2014) 20:2359 Page 11 of 11, 2359
33. Ashbaugh HS, Paulaitis ME (2001) Effect of solute size and
solute-water attractive interactions of hydration water structure
around hydrophobic solutes. J Am Chem Soc 123:10721–10728
34. Christen M, Hünenberger P, Bakowies D, Baron R, Bürgi R,
Geerke DP, Heinz TN, Kastenholz MA, Kräutler V, Oostenbrink
C, Peter C, Trzesniak D, van Gunsteren WF (2005) The GROMOS
software for biomolecular simulation: GROMOS05. J Comput
Chem 26:1719–1751
35. van der Spoel, D van Maaren PJ (2006) The origin of layer
structure artifacts in simulations of liquid water. J Chem Theory
Comput 2:1–11
36. Huang DM, Chandler D (2002) The hydrophobic effect and the
influence of solute-solvent interactions. J Phys Chem B 106:2047–
2053
37. Hummer G, Garde S (1998) Cavity expulsion and weak dewetting
of hydrophobic solutes in water. Phys Rev Lett 80:4193–
4196
38. Kinoshita M (2005) Density and orientational structure of water
around a hydrophobic solute: effects due to the solute size. J Mol
Liq 119:47–54
39. Chiu SW, Clark M, Subramaniam S, Jakobsson E (2000) Collective
motion artifacts arising in long-duration molecular dynamics
simulations. J Comput Chem 21:121–131
40. Harvey SC, Tan RKZ, Cheatham TE (1998) The flying ice cube:
velocity rescaling in molecular dynamics leads to violation of
energy equipartition. J Comput Chem 19:726–740
41. Smith AS (2005) The total solute-water correlation function
forbalancepage Lennard-Jones particles. Fizika A 14:187–
194
42. Leach AR (1996) Molecular modelling: principles and applications.
Addison-Wesley Longman, Harlow
43. Morishita T (2000) Fluctuation formulas in molecular-dynamics
simulations with the weak coupling. J Chem Phys 113(8):2976–
2982