University of Groningen
Effect of Polymer Concentration on the Structure and Dynamics of Short Poly(N,N-
dimethylaminoethyl methacrylate) in Aqueous Solution
Mintis, Dimitris G.; Dompe, Marco; Kamperman, Marleen; Mavrantzas, Vlasis G.
Published in:
Journal of Physical Chemistry B
DOI:
10.1021/acs.jpcb.9b08966
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from
it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date:
2020
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Mintis, D. G., Dompe, M., Kamperman, M., & Mavrantzas, V. G. (2020). Effect of Polymer Concentration on
the Structure and Dynamics of Short Poly(N,N-dimethylaminoethyl methacrylate) in Aqueous Solution: A
Combined Experimental and Molecular Dynamics Study.
Journal of Physical Chemistry B
,
124
(1), 240-252.
https://doi.org/10.1021/acs.jpcb.9b08966
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the
author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license.
More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-
amendment.
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately
and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the
number of authors shown on this cover page is limited to 10 maximum.
Effect of Polymer Concentration on the Structure and Dynamics of
Short Poly(N,N‑dimethylaminoethyl methacrylate) in Aqueous
Solution: A Combined Experimental and Molecular Dynamics Study
Dimitris G. Mintis,
†
Marco Dompe
,
‡
Marleen Kamperman,
§
and Vlasis G. Mavrantzas*
,†,∥
†
Department of Chemical Engineering, University of Patras & FORTH-ICE/HT, GR 26504 Patras, Greece
‡
Physical Chemistry and Soft Matter, Wageningen University, Stippeneng 4, 6708 WE Wageningen, The Netherlands
§
Polymer Science, Zernike Institute for Advanced Materials, University of Groningen, Nijenborgh 4, 9747 AG Groningen, The
Netherlands
∥
Particle Technology Laboratory, Department of Mechanical and Process Engineering, ETH Zu
rich, CH-8092 Zu
rich, Switzerland
*
S
Supporting Information
ABSTRACT: A combined experimental and molecular dynamics (MD) study is
performed to investigate the effect of polymer concentration on the zero shear rate
visc osity η
0
of a salt-free aqueous s olution of poly(N,N-dimet hylaminoethyl
methacrylate) (PDMAEMA), a flexible thermoresponsive weak polyelectrolyte
with a bulky 3-methyl-1,1-diphenylpentyl unit as the terminal group. The study is
carried out at room temperature (T = 298 K) with relatively short PDMAEMA
chains (each containing N = 20 monomers or repeat units) at a fixed degree of
ionization (α
+
= 100%). For the MD simulations, a thorough validation of several
molecular mechanics force fields is first undertaken for assessing their capability to
accurately reproduce the experimental observations and established theoretical laws.
The generalized Amber force field in combination with the restrained electrostatic
potential charge fitting method i s eventually adopted. Three characteristic
concentration regimes are considered: the dilute (from 5 to 10 wt %), the semidilute
(from 10 to 20 wt %), and the concentrated (from 20 to 29 wt %); the latter two are
characterized by polymer concentrations c
p
higher than the characteristic overlap concentration c
p
*. The structural behavior of
the PDMAEMA chains in the solution is assessed by calculating the square root of their mean-square radius of gyration ⟨R
g
2
⟩
0.5
,
the square root of the average square chain end-to-end distance ⟨R
ee
2
⟩
0.5
, the ratio ⟨R
ee
2
⟩/⟨R
g
2
⟩, and the persistence length L
p
.It
is observed that at low polymer concentrations, PDMAEMA chains adopt a stiffer and slightly extended conformation because
of excluded-volume effects (a good solvent is considered in this study) and electrostatic repulsions within the polymer chains.
As the polymer concentration increases above 20 wt %, the PDMAEMA chains adopt more flexible conformations, as the
excluded-volume effects seize and the charge repulsion within the polymer chains subsides. The effect of total polymer
concentration on PDMAEMA chain dynamics in the solution is assessed by calculating the orientational relaxation time τ
c
of
the chain, the center-of-mass diffusion coefficient D, and the zero shear rate viscosity η
0
; the latter is also measured
experimentally here and found to be in excellent agreement with the MD predictions.
1. INTRODUCTION
Poly(N,N-dimethylaminoethyl methacrylate) (PDMAEMA) is
a well-known water-soluble synthetic weak polyelectrolyte
whose pH-responsiveness and thermoresponsiveness are
attributed to the tertiary amine groups on the side chains
along its backbone. Its versatile nature allows for conforma-
tional and dynamical variations upon changing the pH, the
temperature, the ionic strength (salt concentration), the total
polymer concentration, and the molecular weight. For this
reason, it is used in a wide range of applications: as a viscosity
adjusting agent in cosmetics,
1,2
as a solubility enhancer in food
industry,
3,4
and as a drug delivery agent in medical industry.
5
PDMAEMA has been thoroughly studied in the last years
for its potential usage in the formation of complex coacervates
resulting from the complexation of oppositely charged
macromolecules as they undergo an associative liquid−liquid
phase separation.
6
Thermoresponsive complex coacervate
materials play a major role in the development of “high-
performance” underwater adhesi ves because of a unique
combination of properties such as immiscibility with water
and good wetting of the surface.
7
During the last decade,
polyelectrolytes have been greatly utilized in oil and gas
industry as rheology modifiers for controlling and delaying
gelation in the subterranean zone (a reaction zone of
groundwater and seawater) in drilling and fracking processes.
8
Received: September 20, 2019
Revised: December 6, 2019
Published: December 10, 2019
Article
pubs.acs.org/JPCB
Cite This: J. Phys. Chem. B 2020, 124, 240−252
© 2019 American Chemical Society 240 DOI: 10.1021/acs.jpcb.9b08966
J. Phys. Chem. B 2020, 124, 240−252
Downloaded via UNIV GRONINGEN on July 2, 2020 at 17:27:02 (UTC).
See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
Surprisingly, there has been no work, neither experimental
nor computational, addressing the dynamical behavior of
PDMAEMA solutions with a short chain length (N < 50) and
its dependence on polymer concentration. On the contrary,
several previous noteworthy experimental studies
9,10
have
addressed the rheological behavior of PDMAEMA solutions
with considerably longer chain lengths (N > 700); they have
shown that protonated PDMAEMA in solution behaves as a
strong polyelectrolyte but adopts a flexible, coil-like con-
formation in the high-salt (NaCl) limit, both in the semidilute
unentangled and entangled concentration regimes, implying a
shift of its behavior toward that of a neutral polymer. The key
element of the present work is the execution of both
experiments and simulations in order to examine the dynamical
as well as the structural behavior of short PDMAEMA chains
(N = 20) when fully ionized (protonated) in salt-free solutions,
as a function of total polym er concentration at room
temperature (T = 298 K). Our interest in short PDMAEMA
stems from the fact that similar chain length PDMAEMA
chains have been studied experimentally
6,11,12
to determine the
binodal composition of polyelectrolyte complexes of PDMAE-
MA with poly(acrylic acid) (PAA) of practically the same
chain length. Our work here (as well as the recent one on the
molecular dynamics (MD) simulation of aqueous solutions of
short PAA chains
13
) is the first step toward predicting (at a
given pH) the salt−polymer phase diagram for associative
phase separation between PAA a nd PDMAEMA chains
characterized by molecular lengths in the range between 20
and 50 repeat units (i.e., practically identical to those
addressed experimentally) directly from the simulations at
molecular level.
The dimensions of polymers in solutions depend signifi-
cantly on the quality of the solvent and the total polymer
concentration. Graessley
14
has classified polymer solutions into
five different regimes: the dilute regime, the semidilute
unentangled and entangled regimes, and the concentrated
unentangled and entangled regimes. The distinction between
the different regimes has been made based on the variation of
chain dimensions in the solution with total polymer
concentration and molecular weight. Unlike neutral polymers,
the crossover from the dilute to the semidilute concentration
regime for charged polymers (polyelectrolytes) occurs at lower
concentrations. This is primarily attributed to the fact that
charged polymers possess greater chain rigidity in solution than
neutral polymers because of the electrostatic repulsion that
arises between charges, which drastically affects the local
flexibility and tends to increase the global dimensions of the
chain.
15,16
Neutral polymers in a θ solvent retain a coil
conformation,
17
whereas in a good solvent (in the dilute
concentration regime), the coils are expanded because of the
excluded-volume effect.
18
On the contrary, charged macro-
molecules behave as wormlike chains (WLCs), bridging the
polymer behavior from coils to rigid rods.
19,20
For example, in
a recent MD study, PAA, another weak polyelectrolyte, was
found to adopt a WLC conformation when alternately or fully
ionized.
13
Polyelectrolytes are classified as flexible, semiflexible,
and rigid depending on local rigidity (persistence length of the
backbone) and ionic strength (salt concentration).
21,22
The crossover for WLCs from the dilute to the semidilute,
and from the semidilute to the concentrated regime, has been
well defined theoretically by Ying and Chu,
23
by incorporating
two parameters: (a) the effective length L* =(L
p
L
c
)
1/2
, where
L
p
is the persistence length and L
c
is the contour length of the
polymer and (b) an effective diameter d* = d(L
c
/L
p
)
1/4
, where
d is the diameter (thickness) of the macromolecule; in general,
d* is concentration-dependent, as the persistence length is also
concentration-dependent (see also Results and Discussion
section below). The critical overlap concentration of WLCs
from a dilute to a semidilute solution is defined as c
p
* =
(2
3/2
M)/(N
A
L*
3
), where M is the molecular weight of the
polymer and N
A
denotes Avogadro’s number, whereas the
corresponding critical overlap concentration of WLCs from a
semidilute to a concentrated solution is defined as c
p
** =
(0.243M)/(N
A
d*L*
2
). In our study, a comprehensive analysis
of the critical overlap concentration which distinguishes the
dilute from the semidilute and the concentrated regimes is
performed.
The dependence of the chain dimensions of flexible
polyelectrolyte chains on the total polymer concentration has
been described by several scaling theories in the past, as
proposed by De Gennes et al.,
16
Pfeuty,
24
and Dobrynin and
co-workers,
25
in which, however, the polymer persistence
length and other parameters related to chain rigidity (stiffness)
were not considered. In the present study, the global size of
PDMAEMA chains in solution (only a good solvent is
considered) is m easured by estimating the square root
⟨R
g
2
⟩
0.5
of the mean-squared radius of gyration and the square
root ⟨R
ee
2
⟩
0.5
of the mean-squared chain end-to-end distance,
which for simplicity will be denoted in the following as R
g
and
R
ee
, respectively. Dobrynin and co-workers
25
suggest that for
salt-free polyelectrolytes in a good solvent, the value of R
ee
in
dilute solution has no dependence on the total polymer
concentration; thus, R
ee
≈ c
p
0
; however, for c
p
> c
p
*, a scaling
behavior of the form R
ee
≈ c
p
−1/4
is found. The obtained
exponent of −1/4 is consistent with the study of De Gennes
and co-workers
16
andshouldbecontrastedwiththe
corresponding exponent of −1/8 for neutral polymers.
17
A
noteworthy study by Nierlich and co-workers
26
on the mean
radius of gyration of a flexible polyelectrolyte with N = 122 in a
salt-free solution by small-angle neutron scattering revealed
that at a high polymer concentration R
g
decreases as c
p
−1/4
,
whereas at a low polymer concentration, the ratio R
g
/R
g,rod
(with R
g,rod
being the radius of gyration of the macromolecule
if considered to be fully stretched) is around 0.87 and
decreases as the polymer concentration increases. In addition,
an MD study of a linear flexible polyelectrolyte chain modeled
as a freely joined, bead-chain polymer with N = 16, 32, and 64
by Stevens and Kremer
27
revealed that the persistence length
L
p
decreases as the polymer concentration increases, thus
suggesting that at low polymer concentrations the polyelec-
trolyte chain is stiffer than at high polymer concentrations.
Stevens and Kremer
27
also reported that the ratio ⟨R
ee
2
⟩/⟨R
g
2
⟩
at low polymer concentrations is between 8 and 10, implying a
stiff coil conformation, but when the polymer concentration
was increased, they found that ⟨R
ee
2
⟩/⟨R
g
2
⟩ ≅ 6, indicating a
change in chain conformation from being stretched to coiled
(ideal chain).
The effect of polymer concentration on chain dimensions
drastically affects the chain dynamics in the solution. Well-
established theories have been proposed to describe this
dynamics as a function of polymer concentration, such as the
Zimm model for dilute solutions, the Rouse model for
unentangled concentrated solutions, and the Doi−Edwards
model for entangled concentrated solutions and melts.
17,28,29
Dobrynin et al.
25
and Muthukumar
30
found that the dynamics
of linear flexible polyelectrolytes in the dilute concentration
The Journal of Physical Chemistry B Article
DOI: 10.1021/acs.jpcb.9b08966
J. Phys. Chem. B 2020, 124, 240−252
241
regime can be described by the Zimm model because of strong
hydrodynamic interactions. The Zimm model has been found
to be applicable in the case of the dynamics of a weak
polyelectrolyte chain, PAA, at an infinite dilution based both
on experimental
31
and MD
13
studies. According to Dobrynin
and co-workers,
25
in the low-salt limit in the semidilute
unentangled regime, the chain center-of-mass diffusion
coefficient D is concentration-independent, whereas the zero
shear rate viscosity η
0
grows with the square root of the
concentration (η
0
≈ c
p
1/2
), which agrees with the phenom-
enological law proposed by Fuoss.
32,33
In the high-salt
semidilute unentangled regime, the dependence of the zero
shear rate viscosity on polymer concentration is described as η
0
≈ c
p
5/4
. In the case of the low-salt limit in the semidilute
entangled regime, the zero shear rate viscosity scales as η
0
≈
c
p
3/2
, whereas in the high-salt regime, above a certain
concentration where the electrostatic blobs begin to overlap,
the zero shear rate viscosity scales as η
0
≈ c
p
15/4
(i.e., the same
as for uncharged polymers in a good solvent
25
). The diffusion
coefficient D in the low-salt semidilute entangled regime scales
with the total polymer concentration as D ≈ c
p
−1/2
and as D ≈
c
p
−7/4
in the high-salt regime. As described by the theoretical
work of Dobrynin and co-workers,
25
the dependence of
relaxation time τ
c
on polymer concentration for polyelectrolyte
solutions is predicted to decrease (it scales as c
p
−1/2
) in the
unentangled semidilute regime but to be independent of the
polymer concentration in the entangled semidilute regime.
There is still a lack of published theories on the dependency of
the relaxation time on polyelectrolyte concentration in a salt-
free solution in the dilute c oncentration regime. The
experimental work of Wyatt and Liberatore
34
on the solution
rheology of a flexible polyelectrolyte chain (xanthan gum with
a molecular weight of 212 Da) revealed that the relaxation time
τ
c
in the dilute limit scales w ith the total polymer
concentration as τ
c
≈ c
p
1.5
, whereas in the limit above the
characteristic concentration where electrostatic blobs begin to
overlap as τ
c
≈ c
p
4
.
To the best of our knowledge, only one MD simulation
study
35
has addressed the behavior of PDMAEMA chains in
solution. It employed the polymer consistent force field
(PCFF),
36
but no comparison to experiment or theory was
reported. In addition, relatively short simulation times were
accessed, on the order of ∼30 ns. In the present work, we carry
out a systematic analysis of several molecular mechanics force
fields for the behavior of PDMAEMA in aqueous solution and
examine in detail their validity and reliability to reproduce
experimental observations and theoretical laws. Our study
includes a thorough analysis of the critical overlap concen-
tration c
p
* marking the crossover from the dilute to the
semidilute concentration regime for PDMAEMA solutions
with a fixed chain length of N =20atafixed ionization state
(α
+
= 100%) and temperature (T = 298 K). We address the
scaling of chain size and flexibility with the polymer
concentration, including predictions for the diffusion coef-
ficient, the characteristic relaxation time, and the zero shear
rate viscosity, and how they compare to the established theory
and previous studies.
Being a thermoresponsive polymer, the effect of temperature
on the chain dimensions of PDMAEMA would, of course, be
of great importance. However, our main interest in this work is
in the PDMAEMA properties at room temperature because it
is only at this temperature that the equilibrium binodal
composition between the polyelectrolyte complexes or
coacervates of PDMAEMA and fluorescently labeled PAA
with their coexisting dilute phases was measured by Spruijt et
al.,
6,11,12
as a function of polymer chain length and salt
concentration. However, we hope to carry out a systematic
study of temperature effects on the conformation of short
PDMAEMA chains in aqueous solution in the near future, in
order to check how accurately the force field chosen here can
describe the pH-dependence of the lower critical solution
temperature (or even the upper critical solution temperature in
the presence of multivalent counterions in aqueous solution)
in comparison to the available experimental data in the
literature.
37−41
The rest of our work is organized as follows: Section 2
provides a brief description of the experiments carried out and
of the systems simulated and the properties addressed. Section
3 presents a detailed discussion regarding force field validation
and the effect of total polymer concentration on the structural
and dynamical behavior of PDMAEMA in solution. Our paper
concludes with Section 4 summarizing the most important
findings of the present study and briefly highlighting possible
future directions.
Table 1. List of All Simulated Systems Conducted in This Work and Technical Details
no N (mer) c
p
(wt %) c
salt
(M) T (K) force field charge method H
2
O model α
+
(%) total simulation time (ns)
1 30 0 303 MMFF MMFF TIP4P/2005 50 70
2 30 0 338 MMFF MMFF TIP4P/2005 50 70
3 30 0 303 GAFF RESP TIP4P/2005 50 200
4 30 0 338 GAFF RESP TIP4P/2005 50 200
5 30 0 303 OPLS RESP TIP4P/2005 50 200
6 30 0 338 OPLS RESP TIP4P/2005 50 200
7 30 0 303 MMFF RESP TIP4P/2005 50 200
8 30 0 338 MMFF RESP TIP4P/2005 50 200
9 30 0 338 PCFF PCFF PCFF 50 120
10 20 5 0 298 GAFF RESP TIP4P/2005 100 1069
11 20 7 0 298 GAFF RESP TIP4P/2005 100 928
12 20 10 0 298 GAFF RESP TIP4P/2005 100 921
13 20 15 0 298 GAFF RESP TIP4P/2005 100 880
14 20 20 0 298 GAFF RESP TIP4P/2005 100 890
15 20 25 0 298 GAFF RESP TIP4P/2005 100 1002
16 20 29 0 298 GAFF RESP TIP4P/2005 100 1437
The Journal of Physical Chemistry B Article
DOI: 10.1021/acs.jpcb.9b08966
J. Phys. Chem. B 2020, 124, 240−252
242
2. MATERIALS AND METHODS
2.1. Experiments. PDMAEMA with an atactic stereo-
chemistry was purchased from Polymer Source, Inc. Its
number average molecular weight M
n
is 2.7 kg mol
−1
(equivalent to a PDMAEMA chain with the degree of
polymerization N = 20) and its polydispersity index is 1.16.
The polymer was dissolved in Millipore water in concen-
trations spanning from 0 to 29 wt %. The pH of the solution
was adjusted to 5.0.
Rheological measurements were performed on an Anton
Paar MCR501 stress-controlled rheometer equipped with a
Couette geometry (CC-17) with an outer diameter of 18.08
mm, a length of 25 mm, and a bob size of 16.66 mm. The
solution was loaded in the rheometer at a temperature of 20 °C
and then tetradecane was added above the solution to prevent
water evaporation. The sample viscosity η was measured by
performing rotation experiments at shear rates γ
ranging from
100 to 0.01 s
−1
. The torque T decreased as a function of the
shear rate, and the lowest measurable value was 0.01 μN·m.
For every concentration studied, the zero shear rate viscosity η
0
was obtained by averaging the values obtained at the lowest
shear rates, where the Newtonian plateau was obtained, before
reaching the minimum torque.
2.2. MD Simulations. Table 1 provides a list of all systems
simulated in this work by MD together with the detailed
information regarding technical aspects such as the molecular
mechanics force field employed for PDMAEMA chains and
water molecules, charge fitting method used, molecular length
or number of monomers N of the simulated PDMAEMA
chains, simulation temperature T (K), total polyelectrolyte
concentration c
p
(wt %), salt concentration c
salt
(M), degree of
PDMAEMA chain protonation α
+
(%), and total simulation
time (ns).
The chemical structure of the protonated PDMAEMA
molecules assu med in the MD study (relevant to the
experimental chemical structure as reported by Polymer
Source, Inc.) is presented in Figure 1. In all cases, atactic
stereochemistry was adopted.
In addition, and for the purpose of validating the accuracy of
several molecular mechanics force fields on the basis of their
comparison to the predictions made from a previous MD
study,
35
the chemical structure (Figure 2) for the alternate
protonated PDMAEMA molecule was adopted in order to be
consistent with the chemical structure adopted in the previous
study.
The MD simulations were performed using the GROMACS
simulation package version 2016.3.
42−44
All initial con fig-
urations were built in the Materials and Processes Simulations
(MAPS platform, version 4.2, Scienomics SARL, Paris, France)
platform of Scienomics using a modified configurational bias
Monte Carlo scheme.
45,46
A cubic simulation cell with periodic
boundary conditions applied in all directions was used for all
systems. Sufficiently large cubic cells were constructed whose
edge was larger than 2 times the mean end-to-end distance of
PDMAEMA chains in order to prevent spurious interactions of
the molecule with its periodic images in the neighboring cells.
The same energy minimization and equilibration protocol was
applied for all simulated systems. The initial configurations
were energy-minimized by employing the steepest descent
method to eliminate atomic overlaps, with the criterion for
energy convergence set to 50 kJ mol
−1
nm
−1
. Next, a short
simulation of 1 ns in the canonical (nVT) statistical ensemble
was employed (n denotes the total number of interacting units
in the simulation cell, V the volume of the cell, and T the
temperature) with a time step of 1 fs for the integration of the
equations of motion at the temperature of interest (see Table
1). The temperature was maintained constant by using a
Nose
−Hoover thermostat with a coupling constant of 2.5 ps.
Long MD production runs were next carried out in the
isobaric−isothermal (nPT) statistical ensemble where P
denotes the applied pressure (=1 bar), with a time step of 1
fs for the integration of the equations of motion. In this case,
the Nose
−Hoover thermostat was combined with the
Parrinello−Rahman barostat for the simultaneous control of
temperature and pressure, with a coupling constant of 1 ps for
both of them (see Table 1). In the simulations, the bond
lengths and bond bending angles were assumed flexible. The
particle mesh Ewald method
47
was used for computing the
long-range electrostatic interactions.
The pH in our study was adjusted based on the protonation
level of the PDMAEMA molecule. At acidic pH, the
PDMAEMA molecule was assumed to be fully protonated
(all nitrogen atoms were fully protonated); at neutral pH, on
the other hand, the PDMAEMA molecule was assumed to be
alternate protonated (alternate nitrogen atoms were proto-
nated). The case of basic pH was not considered in our study.
The TIP4P/2005 water model was used to solvate the
PDMAEMA chains. Appropriate counterions (Cl
−
) w ere
added to ensure the electroneutrality of the system. Molecule
topologies needed to describe the inter- and intramolecular
interactions of the PDMAEMA chains were created based on
four different force fields: (a) the generalized Amber force field
(GAFF),
48
(b) the optimized potential for liquid simulations
force field (OPLS),
49
(c) the Merck molecular force field
(MMFF), and (d) the PCFF.
36
The GAFF and OPLS
molecule topologies were obtained from the open source
code ACPYPE,
50
the molecule topology for MMFF from the
SwissParam
51
topology generation tool, and the molecule
topology for PCFF from MAPS. The partial atomic charges
were generated based on the restrained electrostatic potential
(RESP)
52
method, with the input for the electronic densities
provided from single-point energy calculations at the Hartree−
Fock level of theory with the 6-31G* basis set. The quantum-
Figure 1. Chemical structure of the protonated PDMAEMA
molecule.
Figure 2. Chemical structure of the alternate protonated PDMAEMA
molecule.
The Journal of Physical Chemistry B Article
DOI: 10.1021/acs.jpcb.9b08966
J. Phys. Chem. B 2020, 124, 240−252
243
