Supporting Information for “Liquid-like states in
micelle-forming diblock copolymer melts”
Kevin D. Dorfman
∗
,
†
and Zhen-Gang Wang
∗
,
‡
Department of Chemical Engineering and Materials Science, University of Minnesota,
Twin–Cities, 421 Washington Ave SE, Minneapolis, MN 55455, USA, and Division of
Chemistry and Chemical Engineering, California Institute of Technology, Pasadena,
California 91125, USA
E-mail: dorfman@umn.edu; zgw@caltech.edu
Contents
1 Detailed Methods
S-3
1.1 Molecular dynamics (MD) simulations . . . . . . . . . . . . . . . . . . . . . S-3
1.2 Conversion into initial guesses for self-consistent field theory (SCFT) . . . . S-4
1.3 SCFT calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S-8
1.4 Structure function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S-11
2 SCFT convergence
S-14
2.1 Body-centered cubic solutions . . . . . . . . . . . . . . . . . . . . . . . . . . S-14
2.2 Liquid-like state solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . S-17
2.3 Computational cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S-19
∗
To whom correspondence should be addressed
†
University of Minnesota, Twin–Cities
‡
California Institute of Technology
S-1
3 Comparison with Frank–Kasper phases
S-22
4 Comparison with the disordered-state free energy
S-23
5 Comparison of bidisperse particle sizes with random particle sizes
S-24
6 Data for different values of
f
A
and
χN
S-28
S-2
1 Detailed Methods
1.1 Molecular dynamics (MD) simulations
The initial conditions for self-consistent field theory (SCFT) calculations were generated by
molecular dynamics (MD) simulations of the Kob–Anderson (KA) fluid
1
using HOOMD–
Blue.
2
The Kob–Anderson fluid is a binary Lennard-Jones fluid consisting of 80% A particles
and 20% B particles, with the particle–particle interaction given by
U
ij
(
r
)
k
B
T
= 4
ij
[
(
σ
ij
r
)
6
−
(
σ
ij
r
)
12
]
(S-1)
The prefactor
ij
and the length scale
σ
ij
are provided in Table S-1. This choice of Lennard-
Jones parameters produces fluids that resist crystallization, leading to a disordered liquid-like
state at higher temperatures and a glass as the temperature is lowered.
1
All MD simulations were performed in the NVT ensemble in a cubic box with periodic
boundary conditions, using a Nos ́e–Hoover thermostat to maintain a dimensionless temper-
ature of 1.5, which produces a liquid-like state. The time step was 0.005 in Lennard-Jones
units. The system was initialized by first placing 1000 particles on a dilute cubic lattice in
a cubic box, disordering them by a short MD simulation, and then slowly compressing the
box to the number density
ρ
MD
= 1
.
20397 of the Kob–Anderson fluid.
1
The dense system
was then equilibrated for 10
6
time steps, followed by a production run of 4
×
10
7
time steps
where the system configuration was sampled every 10
4
time steps.
To generate initial guesses for SCFT from the 4000 frames output from the MD sim-
ulations, we followed the protocol illustrated in Fig. S-1. For a given simulation of 1000
Table S-1: Lennard-Jones parameters in Eq. S-1 for the Kob–Anderson fluid.
Interaction
ij
σ
ij
AA
1.0 1.00
AB
1.5 0.80
BB
0.5 0.88
S-3
Figure S-1: Schematic illustration of the generation of an SCFT initial condition from an
MD simulation frame. The larger A particles are red and the smaller B particles are blue.
Image rendered using fresnel.
3
particles (800 red A particles, 200 blue B particles), we extracted a box from the center of
the simulation volume that would contain 54 particles if the fluid were at the average density
ρ
MD
of the KA fluid. We use separate frames spaced by 10
4
time steps, rather than extract-
ing multiple boxes from the same frame, to ensure that the configurations are uncorrelated;
the computational cost of the MD simulations is small relative to the tens of thousands of
ensuing SCFT calculations. Such a system would thus have the same number of particles as
a 3
×
3
×
3 body-centered cubic (bcc) crystal, which we selected as our large-cell system size
for comparison between the liquid-like states and bcc. Since the KA fluid is compressible,
the number of particles within the central region of the simulation fluctuates, with the sam-
pling producing an approximately Gaussian distribution of particle sizes centered around the
expected mean value of 54 (Fig. S-2a). As expected, samples with a higher number density
tend to also have more small (B) particles (Fig. S-2b).
1.2 Conversion into initial guesses for self-consistent field theory
(SCFT)
The initialization approach described above produced 4000 frames, which were used to gen-
erate independent initial conditions for the SCFT calculations via a modification of the
form-factor approach.
4
Briefly, the volume fraction field for the minority A-block was ini-
S-4
44
46
48
50
52
54
56
58
60
62
64
66
68
Number of particles in KA fluid box
0
.
00
0
.
02
0
.
04
0
.
06
0
.
08
0
.
10
0
.
12
0
.
14
Probability density
44
46
48
50
52
54
56
58
60
62
64
66
68
Number of particles in KA fluid box
0
.
0
2
.
5
5
.
0
7
.
5
10
.
0
12
.
5
15
.
0
17
.
5
Number of small particles
a
b
Figure S-2: Statistics for the number of particles obtained using the approach in Fig. S-1.
(a) Probability density for the number of particles in the central box extracted from the
MD simulation. The snapshots show examples for cases with 46, 54 and 60 particles in the
box. (b) Mean and standard deviation for the number of small (B) particles in the box as a
function of the total number of particles. For 44, 65 and 66 particles, only one sample was
obtained so there is no standard deviation.
tialized in Fourier space by summing over the coordinates
r
j
of each particle
j
,
4
φ
A
(
q
) =
1
V
N
particles
∑
j
=1
f
j
(
q
) exp(
i
q
·
r
j
)
(S-2)
where the form factor for a wavevector of magnitude
q
=
|
q
|
for a sphere of radius
R
j
is
f
j
(
q
) =
4
πR
3
j
[sin(
qR
j
)
−
qR
cos(
qR
j
)]
(
qR
j
)
3
(S-3)
which is formulated with
q
and
R
j
as dimensionless variables in terms of the cubic lattice
vector. In the latter, there is an implicit factor of 2
π
; e.g., for a wavenumber with Miller
indices (
h,k,l
), we define
|
q
|
= 2
π
√
h
2
+
k
2
+
l
2
. To determine
R
j
, the sphere radii were
initialized with a radius
R
j
= 1 to A particles from the KA fluid and
R
j
= 0
.
88 to B particles
from the KA fluid. These radii were then rescaled so that the total volume enclosed by the
spheres corresponds to the total volume fraction of the A-block,
f
A
.
In contrast to prior work,
4
we did not use a Gaussian smearing factor in conjunction with
S-5
a
b
c
Figure S-3: Example of the initialization procedure and refinement during SCFT for
f
A
=
0
.
25 and
χN
= 17.0. (a) Result of the form-factor initialization procedure. (b) Refinement
of the structure after converging SCFT to a tolerance of 10
−
3
. (c) Structure after converging
SCFT to a tolerance of 1
×
10
−
5
. Note the different scales for the color bars as the system
reaches the self-consistent solution.
the form factor to smooth the particle interfaces in the initial guess. The reasons for this
choice are two-fold. First, our calculations tend to be at low values of the segregation strength
χN
, where
χ
is the Flory–Huggins parameter and
N
is the total degree of polymerization.
The resulting interfaces are diffuse, and we did not have any issues converging bcc solutions
from a hard-sphere initial guess. Second, for KA structures, the compressibility of the system
typically leads to overlapping particles when converted to SCFT initial conditions via Eq. S-
2. When applying a smearing factor, we found that the SCFT calculations tended not to
converge when compared to a hard-sphere initial condition, which we attribute to a loss of
the fluid-like structure into a network-like morphology created by the smearing.
To convert the initial guess for
φ
A
into an initial guess for the SCFT fields, we use only
the Flory–Huggins contribution to the chemical potential SCFT fields,
4
ω
A
(
q
) =
χ
(1
−
φ
A
)
ω
B
(
q
) =
χφ
A
(S-4)
initializing the pressure field (Lagrange multiplier) to zero.
The initial guesses for
φ
A
from this protocol for a KA fluid are generally far from the self-
S-6
consistent solution in SCFT, in particular with respect to the incompressibility condition.
Figure S-3a provides one such example. Owing to the compressibility of the KA fluid, there
are many overlapping spheres when we apply the form-factor method, leading to initial
volume fractions that strongly violate incompressibility; note the magnitude of the scale
bar in Fig. S-3a. However, the initial portion of the SCFT calculation (Fig. S-3b) largely
resolves the issues with incompressibility, and the converged solution (Fig. S-3c) has the
desired morphology of a disordered, incompressible liquid.
Unless otherwise noted, calculations for the liquid-like state used this method of ini-
tialization with two particle sizes in the initial condition. The one exception was a set of
calculations discussed in Section 5 at
f
A
= 0.25 and
χN
= 17.0 where we replaced the bidis-
perse sphere sizes in the KA fluid with random sphere sizes to control for artifacts that might
have been introduced by the bidisperse KA model. For the latter, we used the same particle
coordinates obtained from the MD simulations of the KA fluid but assigned each particle a
random radius on the interval [3
/
4
,
4
/
3] rather than the sizes 1.0 and 0.88. The remainder
of the initialization procedure was identical to that described above, including rescaling the
particle radii so that the total minority block volume within the unit cell is
f
A
. Calcu-
lations initialized with this uniform size distribution produce results that are statistically
indistinguishable from those using the bidisperse KA model (see Section 5).
We also attempted to produce solutions using random positions for the particles. These
calculations did not converge, which we attribute to the significant particle overlap generated
by random initialization in a dense state. In principle, initial conditions could also be
generated by simulating a hard-sphere colloid model in a manner similar to what we did here
for the Kob-Anderson fluid, albeit with the somewhat higher computational cost required to
implement the hard sphere model.
S-7
1.3 SCFT calculations
The SCFT equations for a diblock copolymer melt, using the SCFT formulation of Ref. 4,
were solved with the open-source pscf software developed by Morse and coworkers.
4–6
Unless
otherwise noted, the calculations used version 1.0 of the pscf code.
Briefly, we used this code to self-consistently solve the set of polymer field theory equa-
tions for a neat, incompressible diblock copolymer melt with statistical segment length
b
and
total degree of polymerization
N
, where
f
A
=
N
A
/N
is the volume fraction of the minority
A-block and
χN
is the Flory–Huggins parameter. The input to the code are the chemical
potential fields
ω
A
(
r
) and
ω
B
(
r
) for the two species, which were initialized via Eq S-4. The
code solves the modified diffusion equations for the forward propagator
q
(
r
,s
) and reverse
propagator
q
†
(
r
,s
) of a Gaussian chain with a step size
ds
= 0
.
01
N
using the method of
Ranjan
et al.
,
7
and then uses that information to update the fields following the approach
described by Arora
et al.
8
If a self-consistent solution is obtained, the Helmholtz free energy
per chain,
F/nk
B
T
, the optimal values of the lattice parameters in units of the unperturbed
end-to-end distance
R
e
=
bN
1
/
2
, and the volume fraction fields
φ
A
(
r
) and
φ
B
(
r
) are output
by the program. For convenience in later computing the structure factor, we stored the
volume fraction fields in Fourier space, rather than real space.
Two different solution approaches were pursued. In the first, a set of body-centered
cubic (bcc), face-centered cubic (fcc), A15 and
σ
, solutions were computed in a single-unit
cell using symmetry-constrained SCFT with grids of 36
×
36
×
36 for the first three structures
and 128
×
128
×
64 for
σ
. For these calculations, the tolerance based on the definition of
Ref. 8 was set at 10
−
5
. Most of the single unit cell calculations used an earlier version of
pscf and leveraged previously published solutions
9
to generate initial guesses, using unit cells
(cubic or tetragonal) appropriate for the crystal structure. SCFT solutions for the 3
×
3
×
3
(large cell) bcc system and the liquid-like states, unless otherwise noted, were computed
without any imposed symmetry constraints (i.e., using the P1 space group) in a flexible,
orthorhombic unit cell using a 36
×
36
×
36 grid. The latter choice represents a compromise
S-8
45
50
55
60
65
Number of particles
160
180
200
220
Unit cell volume
Figure S-4: Equilibrated unit cell volume, expressed in terms of the cube of the unperturbed
end-to-end distance,
R
3
e
, as a function of the number of particles from the KA fluid initial-
ization for
f
A
= 0.25 and
χN
= 17.0. For the same conditions, the large-cell bcc calculation
produces a unit cell volume of 169.15
R
3
e
.
between the speed of the calculation, memory required to store the solutions, and the desired
resolution for the free energy, which we will address in Section 2. The large-cell calculations
were converged to a tolerance of 10
−
5
, consistent with that used in prior work for Frank–
Kasper phases,
9
which have similarly large numbers of particles per unit cell (but a higher
symmetry). To achieve this value of the tolerance, we needed to increase the penalty for unit
cell stress from its standard multiplier of 10 to a value of 50. This change is implemented
using the
scaleStress
option in pscf.
6
All large-cell calculations used a flexible orthorhombic box whose volume was adjusted
in conjunction with the solution of the SCFT equations to minimize the free energy.
8,10
This approach allows the system to exchange mass and adjust the distribution of micelle
aggregation numbers; this change is evident in Fig. S-3 in the maximum value of the scale
bars, which represents the volume fraction of the minority block in the centers of the par-
ticles. The typical asymmetry in the lattice parameters was a few percent, which reflects
the compromise between the isotropic nature of the macroscopic liquid and the finite size
of the simulation box. As anticipated, the equilibrium unit cell volume tends to increase
with increasing number of particles from the KA initialization (Fig. S-4), although there
S-9
is considerable overlap between unit cell volumes even when different numbers of particles
are used for initialization. The general trend in Fig. S-4 holds for all values of
f
A
and
χN
studied here (see Figs. S-19 to S-22).
The SCFT equations were solved using Anderson mixing.
11
While the large-cell bcc cal-
culations converged easily, the convergence starting from initial guesses produced by the KA
fluid proved to be challenging for three reasons. The first issue is already apparent from
Fig. S-3, i.e., that the initial guesses tend to be relatively far from the saddle point, in par-
ticular for the incompressibility condition that corresponds to the saddle maximum on the
imaginary axis. Thus, a non-trivial number of iterations are required to move the system
closer to incompressibility, as illustrated by the scale bars in Fig. S-3. The second issue
arises from the lack of symmetry, which leads to all of the Fourier modes being important
for the solution to the problem. The resulting optimization problem is more challenging than
the large-cell bcc problem, where the noise from non-zero symmetry-forbidden wavenumbers
is ultimately attenuated during the iterative process, and much more challenging than a
symmetry-constrained solution for bcc, where symmetry relationships lead to a large reduc-
tion in the number of unknowns for a given grid size. The third issue was less obvious: often,
we found that simply trying to converge a KA fluid initial guess to an arbitrary tolerance
value can lead to a stalled iterator. Interestingly, restarting the calculation (i.e., deleting the
Anderson mixing history) often allowed the solution to escape from the stall. Both the sec-
ond and third issues suggest challenges with Anderson mixing for this problem; investigating
this numerical issue is a promising avenue of future research.
To circumvent these issues within the pscf software, the large-cell calculations were per-
formed by setting 200 logarithmically-spaced target values for the convergence between an
upper bound of 10
−
2
and a lower bound 10
−
5
; the small gap between tolerance values as we
get closer to the lower bound facilitated slowly moving solutions to reach the target toler-
ance. Moreover, there are points in the solution where the Anderson mixing iterator can
drift between larger and smaller values of the tolerance, which would ultimately lead to a
S-10
stalled solution. The logarithmically-spaced target tolerances act like “adsorbing boundary
conditions” for the iterator, allowing it to capture a solution at a lower tolerance value and
then restart the calculation. In addition to helping alleviate the problem of stalled solutions,
this step-wise approach outputs values of the free energy as a function of the tolerance used
for that calculation, which allowed us to assess convergence of the free energy in addition
to convergence for self-consistency without any additional modifications to the pscf source
code.
1.4 Structure function
Calculation of the structure function
S
(
q
) was facilitated by the
kgrid
output of the pscf
v1.0 software, which provides the output of FFTW package for the volume fractions
φ
A
and
φ
B
in terms of the Miller indices (
h,k,l
). To compute the structure function, we first
computed
S
(
q
) =
∑
q
=2
π
√
h
2
+
k
2
+
l
2
q
6
=0
φ
A
(
h,k,l
)
φ
A
(
−
h,
−
k,
−
l
)
(S-5)
where the index for the summation corresponds to all Miller indices with common values
of
q
. This dimensionless result was converted to physical units by noting that the SCFT
problem is formulated in terms of the unperturbed end-to-end distance
R
e
=
bN
1
/
2
. We
thus convert this initial calculation of
S
(
q
) into the experimentally-relevant form
S
(
qR
g
) by
dividing by rescaling (
h,k,l
) by their lattice parameter and then dividing by a factor of
√
6,
the latter accounting for the difference between
R
e
and the radius of gyration,
R
g
.
For bcc, this approach furnishes the Bragg peaks in the structure factor. Figure S-5
provides an example of a bcc solution for
f
A
= 0.25 and
χN
= 17.0, illustrating both the
allowed reflections for bcc and the impact of the micellar form factor on the peak intensities.
For the liquid-like states, the initial analysis produces Bragg peaks as well (e.g., Fig. S-6,
reproduced from the main text here for convenience) due to the periodic boundary condition
on the calculation cell. Owing to the flexible orthorhombic unit cell, each grid point con-
S-11
0
2
4
6
8
qR
g
10
−
7
10
−
6
10
−
5
10
−
4
10
−
3
10
−
2
10
−
1
S
(
qR
g
)
√
1
√
2
√
3
√
4
√
5
√
6
√
7
√
8
√
9
√
10
√
11
√
12
√
19
Figure S-5: Structure factor for bcc for
f
A
= 0.25 and
χN
= 17.0. The inset shows the
morphology of the SCFT solution.
0
2
4
6
qR
g
10
−
8
10
−
7
10
−
6
10
−
5
10
−
4
10
−
3
S
(
qR
g
)
Figure S-6: Structure factor for one KA initial condition for
f
A
= 0.25 and
χN
= 17.0. This
plot also appears in the main text.
tributes a Bragg peak to the structure factor because the different lattice parameters (
a,b,c
)
break the degeneracy in (
h,k,l
) that would otherwise be present in a cubic unit cell, leading
to the dense set of peaks in Fig. S-6. Moreover, each initial condition provides a different
set of Bragg peaks, owing to both the packing of the particles within that system and the
different unit cell lattice parameters.
We thus treated each of the scattering patterns obtained for a given converged KA initial
guesses as representative of one grain within a macroscopic sample. In such a system, we
then anticipate that
S
(
q
) can be approximated by the average over all of the Bragg-like
S
(
q
)
S-12
0
2
4
6
qR
g
10
−
6
10
−
5
10
−
4
10
−
3
10
−
2
10
−
1
S
(
qR
g
)
Figure S-7: Structure function obtained for
f
A
= 0.25 and
χN
= 17.0 after applying Gaussian
smearing to each Bragg peak and then averaging over all of the converged solutions.
patterns. To assist with the averaging, we applied a Gaussian smoothing to the discrete
patterns like Fig. S-6 by converting peak
k
at
q
k
to a smeared value
̃
S
(
q
;
q
k
) =
S
(
q
k
)
(2
πσ
2
)
1
/
2
exp
[
−
(
q
−
q
k
)
2
2
σ
2
]
(S-6)
where
q
is a grid of 5000 points between
q
= 0 and a value of
qR
g
that is slightly larger than
the maximum value of
q
with a Bragg peak. To smooth the data, we use
σ
= 0.01, summed
the
̃
S
(
q
;
q
k
) over
q
k
to get a smoothed structure factor for a given converged solution, and
then averaged all of the smoothed structure factors to get the ensemble average. Figure S-7
shows the result of smoothing all of the results from the data set that produced the example
in Fig. S-6. The smoothed scattering patterns for all conditions are provided in Figs. S-23
to S-26.
To avoid a proliferation of symbols, both the Bragg peaks produced by Eq. S-5 and the
Gaussian smoothing of Eq. S-6 and subsequent averaging of all converged solutions are both
denoted as
S
(
qR
g
); it should be clear from the figure whether the data reported as
S
(
qR
g
)
are the discrete Bragg peaks of Eq. S-5 or the average of smeared peaks given by Eq. S-6.
To compute the peak in the structure factor,
q
∗
, we used a block averaging approach
with 64, 32,
...
, 2 blocks. The uncertainty in
q
∗
was estimated as the optimal blocking that
S-13