Supplementary Information: Enhancing Cation
Diffusion and Suppressing Anion Diffusion via
Lewis-Acidic Polymer Electrolytes
Brett M. Savoie, Michael A. Webb, Thomas F. Miller III
Division of Chemistry and Chemical Engineering, California Institute of Technology
Pasadena, California 91125, United States
January 8, 2017
Contents
1 Computational Methods
S2
1.1 Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S2
1.2 Ion Diffusion Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S3
1.3 Solvation Free-Energy Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . S4
1.4 Salt Lattice Energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S6
2 Results Referenced in the Main Text
S6
2.1 Dilute Simulations for All Ions and Polymers . . . . . . . . . . . . . . . . . . . . . . S6
2.2 Finite Salt Concentration Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . S6
3 Parametrization Details
S11
3.1 Force-Field Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S11
3.2 Fitting Procedure for Intramolecular Modes . . . . . . . . . . . . . . . . . . . . . . . S11
3.3 Fitting Procedure for Lennard-Jones Parameters and Atomic Charges . . . . . . . . S12
4 Force-Field Validation for PEO
S15
5 Force-Field Parameter Tables
S17
5.1 Syntax for Defining Atom Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S17
5.2 PEO Parameter Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S18
5.3 CBC Parameter Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S20
5.4 CBCC Parameter Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S22
5.5 HBC Parameter Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S24
5.6 HBCC Parameter Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S26
5.7 Cl
-
Parameter Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S28
5.8 Triflate Parameter Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S29
5.9 TFSI
-
Parameter Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S30
S1
6 Force-Field Fit Potentials: Bonds, Angles, Dihedrals
S31
6.1 PEO Bond, Angle, and Dihedral Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . S31
6.2 CBC Bond, Angle, and Dihedral Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . S33
6.3 CBCC Bond, Angle, and Dihedral Fits . . . . . . . . . . . . . . . . . . . . . . . . . . S36
6.4 HBC Bond, Angle, and Dihedral Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . S38
6.5 HBCC Bond, Angle, and Dihedral Fits . . . . . . . . . . . . . . . . . . . . . . . . . . S41
6.6 Triflate Bond, Angle, and Dihedral Fits . . . . . . . . . . . . . . . . . . . . . . . . . S44
6.7 TFSI
-
Bond, Angle, and Dihedral Fits . . . . . . . . . . . . . . . . . . . . . . . . . . S46
7 Force-Field Fit Potentials: Lennard-Jones Parameters
S49
7.1 PEO Lennard-Jones Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S49
7.2 CBC Lennard-Jones Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S51
7.3 CBCC Lennard-Jones Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S53
7.4 HBC Lennard-Jones Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S55
7.5 HBCC Lennard-Jones Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S57
7.6 Li-Anion Lennard-Jones Fits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S59
1 Computational Methods
1.1 Molecular Dynamics
The molecular dynamics (MD) simulations employed a non-polarizable force-field that was sys-
tematically parameterized using the Topologically Automated Force-Field Interactions (TAFFI)
methodology; a summary of the methodology is provided here with additional comprehensive de-
tails provided in the later sections.
All intramolecular modes, atomic partial charges, and ion-polymer Lennard-Jones parameters
were parameterized in this work to provide a consistent description of all polymers. Intramolecular
modes, such as bonds, angles, and dihedrals, were parameterized by fitting standard harmonic and
fourier potential energy terms to mode-scans performed at the B3LYP-D3/def2-TZVP level. For
each polymer, a model compound comprised of a tetramer of the polymer capped with methyl
groups was employed for the intramolecular mode parametrization. Approximate atomic partial
charges were obtained from CHELPG calculations based on the B3LYP-D3/def2-TZVP electron
densities at the optimized geometries. After obtaining the intramolecular modes and approximate
partial charges, condensed phase MD simulations were performed on the pure oligomers and also
solutions with each salt to provide configurations to further refine the atomic partial charges and
parametrize the dispersion interactions. Final partial charges on each atom were determined by
averaging the results of CHELPG calculations based on the B3LYP-D3/def2-TZVP electron den-
sities for one hundred molecular configurations sampled from condensed phase MD trajectories.
Pairwise molecular configurations sampled from the condensed phase MD trajectories were used
to parametrize the dispersion contributions to the polymer-polymer, ion-polymer, and ion-ion in-
teractions. The dispersion interactions were parameterized by fitting Lennard-Jones potentials to
the residual of the fixed electrostatic interactions and the counter-poise corrected B3LYP-D3/def2-
TZVP interaction energies calculated for the pairwise configurations. Full parametrization details,
parameter tables, mode scans, and fit potentials are included in later sections.
All MD simulations were performed within the LAMMPS software suite.
1
All trajectories em-
ployed periodic boundary conditions, particle-particle-particle-mesh (pppm) evaluations of long-
range interactions beyond a 14
̊
A cut-off, a Nos ́e-Hoover barostat with 1000 fs relaxation, and a
Nos ́e-Hoover thermostat with 100 fs relaxation (NPT). Equations of motion were evolved using the
S2
velocity-Verlet integrator and a two femtosecond timestep for polymers without explicit hydrogen
atoms and a one femtosecond timestep for polymers with explicit hydrogen atoms. Intramolecular
pairwise interactions for atom pairs connected by fewer than four bonds were excluded during the
MD simulations to avoid double counting with dihedral energy terms. Ion diffusion and solvation
free energy trajectories were initialized from a common set of four independently equilibrated neat
polymer trajectories. Each neat polymer trajectory included a single polymer chain with a mass
of approximately 30 kDa that was initialized using a protocol to randomize chain orientation and
avoid configurations with significant overlap between atoms. These configurations were initially
equilibrated at a temperature of 10 K and a pressure of 50 atm for 50 ps, followed by a 10 ns
annealing phase at a temperature of 500 K and a pressure of 1 atm. After annealing, the config-
urations were simulated for 11 ns at a temperature of 400 K and a pressure of 1 atm to collect
production data. The first nanosecond of the production trajectories was used for equilibration and
the remaining 10 ns were used to confirm the convergence of the density. The final configurations
from the neat polymer trajectories were used as input geometries for the ion diffusion and solvation
free energy trajectories.
1.2 Ion Diffusion Simulations
For each combination of ion and polymer, sixteen independent ion-diffusion trajectories were per-
formed. These trajectories were initialized from the four pre-equilibrated neat polymer trajectories
by randomly inserting a single ion into each configuration (four independent insertions for each
pre-equilibrated configuration). Each trajectory included a single ion to avoid correlated ion mo-
tions, and the excess charge was neutralized with a uniformly distributed background charge. The
initial geometry was relaxed by performing 1 ps of NVE dynamics with atom displacements limited
to 0.01
̊
A for each timestep, followed by 1 ns of NPT dynamics at a temperature of 400 K and
pressure of 1 atm. After relaxation, an additional 300 ns of NPT dynamics were performed to
collect production data. A total of 4.8
μ
s of ion diffusion dynamics were collected for each ion in
each polymer (16 independent trajectories, each 300 ns long).
Diffusivities can be calculated from the long-timescale trajectories using the Einstein equation
D
i
= lim
t
→∞
d
〈
|
r
i
(
t
)
−
r
i
(0)
|
2
〉
6d
t
,
(1)
where
D
i
is the diffusion coefficient for ion,
i
, and the term in brackets is the MSD evaluated at
time
t
. Because many of the systems studied still exhibit sub-diffusive behavior even at long times,
apparent diffusivities are reported by approximating the derivative in Eq. (1) by finite difference
using the MSD at
t
= 150 and
t
= 0 ns.
2, 3
The Li
+
transference number,
T
Li
, was calculated using
T
Li
=
D
Li
/
(
D
Li
+
D
anion
).
Contact durations were calculated by defining a characteristic function,
h
ij
(
t
), that reports on
contacts between pairs of atoms
h
ij
(
r
) =
{
1
,
for
r
ij
≤
r
T
0
,
for
r
ij
> r
T
,
(2)
where
i
and
j
denote atomic indices, and
r
T
is a pair-specific threshold distance, which is chosen
based on the size of the first coordination shell for the corresponding atom types. Specifically,
r
T
=
r
max
+ 2
σ
, where
r
max
is the radial separation at the first maximum of the corresponding ion-
polymer radial-distribution function (see Supporting Information), and
σ
is the standard deviation
obtained by fitting the full width at half maximum of the peak to a Gaussian function. The
S3
time autocorrelation function,
〈
h
ij
(0)
h
ij
(
t
)
〉
, was calculated for each contact type by averaging
the trajectories with respect to time. Standard errors were calculated by separately averaging the
results from the sixteen independent trajectories.
1.3 Solvation Free-Energy Calculations
Thermodynamic integration was used to calculate the ion-specific solvation free energies in each
polymer.
4
For each combination of ion and polymer, sixteen independent ion-insertion trajectories
were used. Scaled potentials were used to gradually introduce the ion-polymer potential energy
terms, and convergence was facilitated by introducing the polymer-ion Lennard-Jones interactions
before introducing the polymer-ion electrostatic interactions. A scaled potential was used to first
introduce the polymer-ion Lennard-Jones interactions
U
LJ
(
λ
1
) =
U
P
+
λ
1
(
U
P+X
LJ
−
U
P
)
,
(3)
where
λ
1
is the scaling parameter,
U
P
is the potential of the pure polymer, and
U
P+X
LJ
is the
potential of the pure polymer plus ion-polymer Lennard-Jones interactions. Standard
λ
-dependent
soft-core Lennard-Jones potentials, as implemented in LAMMPS with
n
= 1 and
α
LJ
= 0
.
5, were
used for all ion-polymer interactions to smooth the potential energy function.
5
When
λ
1
= 0, the
ion and polymer are non-interacting, and when
λ
1
= 1, the ion-polymer Lennard-Jones terms fully
contribute to the potential energy. A second scaled potential was used to subsequently introduce
the polymer-ion electrostatic interactions
U
C
(
λ
2
) =
U
P+X
LJ
+
λ
2
(
U
P+X
LJ
+X
C
−
U
P+X
LJ
)
,
(4)
where
U
P+X
LJ
+X
C
is the potential of the polymer plus all polymer-ion interactions. The potential
in Eq. (4) was implemented by using
λ
-scaled charges on the ion. When
λ
2
= 0, the ion and
polymer interact only through the Lennard-Jones potential, and when
λ
2
= 1, both the ion-polymer
electrostatic and Lennard-Jones terms fully contribute to the potential energy. The total solvation
free-energy was obtained by
∆
G
TI
=
∫
1
0
〈
d
U
LJ
d
λ
1
〉
d
λ
1
+
∫
1
0
〈
d
U
C
d
λ
2
〉
d
λ
2
.
(5)
The brackets in Eq. (5) indicate an ensemble average, and the approximation has been made that
the
P
∆
V
contribution to the free energy change can be safely neglected. The integrals in Eq. (5)
were evaluated numerically using the trapezoidal rule, with
λ
1
and
λ
2
incremented in steps of 0.1
(twenty-one steps total, eleven for the Lennard-Jones phase and eleven for the electrostatics, less one
redundant step connecting the two phases). The system was allowed to equilibrate for 100 ps at each
λ
-step, then an additional 100 ps of dynamics were used for calculating the necessary derivatives.
The derivatives in Eq. (5) were calculated by finite-difference. At endpoints, forward or backward
finite-difference was used, at all other points the central difference was used with a
λ
-step of 0.01 to
evaluate the derivative. In the case of TFSI
−
, an additional free energy contribution associated with
removing the intramolecular electrostatics must be computed. Free-energy perturbation was used
to evaluate this contribution from a ten ns MD trajectory of a single TFSI
−
molecule in vacuum.
The reported ∆
G
TI
values were calculated as the average over all ion-insertion trajectories, with
errors in the mean estimated by bootstrap resampling (5 million samples).
6
It is common for borane centers to undergo changes in bonding hybridization and associated
structural rearrangements upon anion complexation.
7
Since these effects are not captured by the
force-field used in this work, it is anticipated that the MD thermodynamic integration calculations
S4
significantly underestimate the anion solvation free energies in the Lewis-acidic polymers. Using
ab
initio
geometry optimizations and energy calculations, we thus include a free-energy-perturbation
correction to the solvation free energy (for both the Lewis-acidic polymers and PEO) associated
with local structural relaxation during anion complexation,
∆
G
corr
= ∆
G
rel
,
anion
−
∆
G
rel
,
neat
,
(6)
where
∆
G
rel
,
anion
=
−
β
−
1
ln
〈
e
−
β
∆
U
rel
,
anion
〉
(7)
accounts for the free-energy change in the anion-polymer system upon relaxation, and
∆
G
rel
,
neat
=
−
β
−
1
ln
〈
e
−
β
∆
U
rel
,
neat
〉
(8)
accounts for the corresponding effect in the neat polymer. The angle brackets in Eqs. (7) and
(8) correspond to ensemble averaging over the MD configurations, and the configuration-specific
relaxation energies ∆
U
rel
,
anion
and ∆
U
rel
,
neat
are obtained from quantum chemistry calculations at
the PBE-D3/def2-SVP level of theory, as described below.
The ∆
U
rel
,
anion
terms in Eq. (7) were calculated by performing constrained geometry optimiza-
tions of the anion solvation structures. For each anion in each polymer, ten solvation snapshots
were extracted at 100-picosecond intervals from the dilute ion concentration MD trajectories. Each
solvation snapshot consisted of the anion and all polymer atoms within
r
outer
= 5
̊
A of any anion
atom. The value of
r
outer
was chosen to be large enough to include the first solvation shell for
all combinations of polymer and anion (
r
inner
= 4
.
5
̊
A, see Fig. S1) plus a buffer region of 0
.
5
̊
A.
The united-atom moieties and the terminal atoms of polymer fragments in the solvation snapshot
were then hydrogenated, and a constrained geometry optimization of the added hydrogens was
performed. Subsequently, a constrained geometry optimization was performed in which all binding
atoms (carbon for PEO and boron for all Lewis-acidic polymers) within
r
inner
of any anion atom
and polymer atoms within two bonds of those binding atoms were relaxed, while the positions of
all anion atoms and remaining polymer atoms were constrained. For instances where these cutoff
values return a polymer fragment without any constrained atoms, a position constraint was placed
on the atom in the fragment farthest from the anion. The relaxation energy, ∆
U
rel
,
anion
, is obtained
from the single-point energy difference for the system before and after the constrained geometry
optimization. For each combination of polymer and anion, the average number of polymer binding
atoms that were within
r
inner
,
N
poly
, is determined, and this value is used in the protocol for the
calculation of ∆
U
rel
,
neat
.
The ∆
U
rel
,
neat
terms in Eq. (8) were similarly calculated by performing constrained geome-
try optimizations of snapshots obtained from neat-polymer MD trajectories. The procedures for
sampling and optimizing the neat polymer snapshots were identical to the solvation geometries,
except that the neat polymer snapshots included all polymer atoms within
r
outer
= 7
.
5
̊
A of the
center of the simulation box; this larger value of
r
outer
was used in the neat-polymer calculations
to account for the reduced density of binding atoms in the neat-polymer snapshots relative to the
anion-containing snapshots. During the geometry relaxation, only the
N
poly
binding atoms closest
to the center of the simulation box and the polymer atoms within two bonds of those binding atoms
were allowed to relax, where
N
poly
is defined above; this ensures that the same average number
of binding atoms are relaxing in calculating both ∆
U
rel
,
anion
and ∆
U
rel
,
neat
. The neat relaxation
energy, ∆
U
rel
,
neat
, is obtained from the single-point energy difference for the system before and
after the constrained geometry optimization.
S5
The final expression used to evaluate the solvation free energy for an anion was
∆
G
S
= ∆
G
TI
+ ∆
G
corr
.
(9)
Table S1 includes the values of ∆
G
corr
for all anions in all polymers. For Li
+
in all polymers,
∆
G
S
= ∆
G
TI
.
1.4 Salt Lattice Energies
The salt lattice energies for LiCl and LiTriflate presented in Fig. 5b were taken from reference 16.
No experimental reference exists for LiTFSI, so the lattice energy was estimated using a modified
Kapustinki approach
8
developed by Jenkins and a volume of 136.1
̊
A
3
for TFSI calculated based
on the DFT optimized structure at the B3LYP-D3/def2-TZVP level.
2 Results Referenced in the Main Text
2.1 Dilute Simulations for All Ions and Polymers
Fig. S1 shows the ion-polymer radial distribution functions
g
AB
(
r
) for all ions in all polymers
simulated in the dilute-ion regime. For each ion-polymer combination,
g
AB
(
r
) is calculated such
that
A
refers to the set of atoms associated with the ion and
B
refers to the set of atoms associated
with the polymer that have partial charges with opposite sign as the ion.
Fig. S2 shows the mean-squared displacement (MSD) for all ions, primary ion-polymer contact
durations, and the MSDs of the polymer segments calculated from the dilute-ion simulations at
400 K. The data for Li
+
and Cl
–
are reproduced from the main text for comparison.
2.2 Finite Salt Concentration Simulations
All finite salt concentration simulations were initialized using the same equilibrated neat polymer
configurations that were used to initialize the dilute-ion simulations, as described in the main text.
Ions were added to the neat polymer configurations at random positions until reaching the specified
concentration. The initial geometry was relaxed by performing 5,000 energy minimization steps,
during which the displacements of each atom were limited to 0.01
̊
A per step, followed by 5 ns of
NPT dynamics at a temperature of 400 K and pressure of 1 atm. After relaxation, NPT dynamics
at a temperature of 400K and pressure of 1 atm were run for an additional 75 ns (using a 1 fs
timestep) for polymers with explicit hydrogen atoms and 150 ns (using a 2 fs timestep) for all other
polymers. A Nos ́e-Hoover thermostat (100 fs relaxation timescale) and barostat (1000 fs relaxation
timescale) were used to control the temperature and pressure for all NPT simulations. Sixteen
independent trajectories were simulated per electrolyte at each concentration.
The ion force-fields were parametrized in the dilute limit and are not expected to be quantita-
tively valid at finite salt concentrations, where polarizability plays an important role. Simulations of
polymer electrolytes at high salt concentrations have also been shown to exhibit increased aggrega-
tion and suppressed conductivity relative to experiment when using non-polarizable or fixed-charge
force fields.
9–13
Likewise, initial simulations with full charges on the ions showed negligible con-
ductivity and failed to reach equilibrium even at long (75 ns) timescales (not shown). To obtain
reduced pairing and mobile ion dynamics, the charges on the ions were uniformly scaled by 0.5
during the simulation, similarly to previous work.
9–12
Any attempt to compare the ion-pairing
results in Fig. S3 with the dilute ion concentration solvation free energies in Fig. 5 of the main
S6
text should consider that the finite-concentration results in Fig. S3 are simulated with scaled point
charges and do not include intramolecular relaxation corrections like those used in Fig. 5.
Fig. S3 shows the conductivity (
σ
), degree of uncorrelated ion motion (
α
), and transference
number (
T
Li
) for each polymer with LiTFSI at two concentrations (
r
=0.1 and
r
=0.05, where
r
is the ratio of Li
+
ions to polymer monomer units) and scaled partial charges on the salt. The
conductivity was calculated according to
σ
= lim
t
→∞
e
2
6
tV k
B
T
N
∑
ij
z
i
z
j
〈
d
i
(
t
)
d
j
(
t
)
〉
,
(10)
where
e
is the fundamental charge,
V
is the simulation volume,
T
is the temperature,
z
is the
integer charge of the ion,
d
is the vectorial displacement of the ion at time
t
, the
i
and
j
indices
run over all ions in the simulation, and
N
is the total number of ions in the simulation box. When
calculating
σ
, the charges on all ions were unscaled.
The degree of uncorrelated ion motion,
α
, was calculated according to
α
= lim
t
→∞
1
6
t
(
D
Li
+
D
anion
)
N
N
∑
ij
z
i
z
j
〈
d
i
(
t
)
d
j
(
t
)
〉
,
(11)
where
D
Li
and
D
anion
correspond to average the diffusivity of Li
+
and the anion, respectively,
calculated as
D
β
= lim
t
→∞
1
6
t
∑
i
∈
β
〈
d
i
(
t
)
d
i
(
t
)
〉
,
(12)
and the summation in Eq. 3 runs over all ions of type
β
. When the cross terms in the summation
in Eq. 2 average to zero (i.e., uncorrelated ion motion), then
α
equals one. When Li
+
motion
perfectly correlates with anion motion (i.e. paired behavior),
α
equals zero. As in the main text,
T
Li
was calculated according to
T
Li
=
D
Li
D
Li
+
D
anion
.
(13)
The limit on
t
was taken as 75 ns when calculating all quantities. The data shown in Fig. S3 repre-
sents averages over all trajectories (16 trajectories per electrolyte), with the error bars representing
the standard error for each quantity calculated across all trajectories.
Table S1: Free energy corrections, ∆
G
corr
, associated with structural relaxations during anion
complexation for all anions in all polymers.
Cl
−
Triflate
TFSI
PEO
−
10
±
4
−
14
±
9
−
1
±
6
CBC
−
37
±
6
−
19
±
6
−
74
±
14
CBCC
−
100
±
9
−
28
±
8
−
26
±
5
HBC
−
83
±
14
−
10
±
10
−
33
±
14
HBCC
−
67
±
13
−
41
±
15
−
61
±
8
Units in kcal
·
mol
−
1
.
Standard errors are reported.
S7
0
2
4
6
8
10
0
5
10
15
20
25
30
CBCC
CBC
PEO
HBCC
HBC
0
2
4
6
8
10
0.0
0.5
1.0
1.5
2.0
0.0
0.5
1.0
1.5
2.0
2.5
Separation (
Å
)
g
Li
+
X
-
(r)
g
Cl
-
X
+
(r)
g
Triflate
-
X
+
(r)
g
TFSI
-
X
+
(r)
Figure S1: Comparison of radial-distribution functions (RDFs) for all ions in all polymers simulated
in the dilute-ion regime at 400 K. “X” in each case includes all polymer atoms of opposite partial
charge as the ion. The data are represented as a histogram with a bin width of 0.1
̊
A for
r
.
S8