Monte Carlo Simulations of Diblock Copolymers with an Interactive Front-End

Diblock Copolymers (DBCs)

A diblock copolymer consists of an A-block of $N_A$ segments and a B-block of $N_B$ segments (figure 1), joined together at one end to give a linear chain with $N=N_A+N_B$ total segments and composition fraction $f_A=N_A/N$.




The strength of interactions between A and B segments is controlled by the Flory-Huggins parameter, $\chi$, which is inversely related to temperature (i.e., low $\chi$ corresponds to high temperature). When the A-B interactions in a melt of DBCs become sufficiently unfavourable (i.e., $\chi$ becomes large enough) they microphase-separate, forming periodically ordered structures with domains on the nanometer scale. The equilibrium structures include the classical lamellar (L), cylindrical (C) and spherical (S) phases, as well as the complex gyroid (G) and Fddd phases (see figure 2).

Self-consistent Field Theory (SCFT)

Imagine looking at a single DBC in a melt. It has complicated interactions with all the other polymers in the system and capturing them all is a difficult mathematical and computational challenge. SCFT addresses this problem by replacing the interactions that a polymer has with other chains by mean fields.

Monodisperse Phase Diagram
SCFT has been enormously successful in the study of DBCs, producing the famous equilibrium phase diagram for a monodisperse melt shown in figure 3. It displays the various ordered phases, which are separated from the disordered phase by the order-disorder transition (ODT). It is also worth noting that the phase diagram is symmetric about $f_A=1/2$.

Limitations of SCFT

Unfortunately, the mean-field nature of SCFT means that it doesn't include the fluctuation effects that are inherent in real statistical systems. As a result, we are interested in fluctuation corrections to figure 3, which are controlled by the invariant polymerisation index, $\overline{N}=a^6{\rho }_0^{2}N$, where ${\rho }_0$ is the segment density and $a$ is the statistical segment length.

SCFT corresponds to the case where $\overline{N}=\mathrm{\infty }$, with fluctuations increasing as $\overline{N}$ is reduced. The main effects of fluctuations on the monodisperse phase diagram in figure 3 are that the ODT gets shifted upward, and we get direct transitions between the disordered and ordered phases. For example, there will be a direct transition beween gyroid and disorder as indicated schematically in figure 4.

Monodisperse Phase Diagram from Experiments

In 1994, Bates et. al. performed experiments to produce monodisperse phase diagrams for three different DBC chemistries: PE-PEP ($\overline{N}=27\times {10}^3$), PEP-PEE ($\overline{N}=3.4\times {10}^3$) and PI-PS ($\overline{N}=1.1\times {10}^3$). Each of the phase diagrams exhibited the spherical, cylindrical and lamellar phases, but no gyroid or Fddd. There was, however, a region of another phase known as perforated lamellar in the region where SCFT predicts gyroid. Indeed, later experiments by Hajduk et. al. and Vigild et. al. showed that perforated lamellar is actually a long-lived metastable state that eventually converts to gyroid. Even more recent experiments by Wang et. al. and Takenaka et. al. have also confirmed the presence of Fddd on both sides of the phase diagram.

Whilst the experiments and SCFT phase diagrams exhibit a lot of good qualitative agreement, the significant differences are due to the effects of fluctuations near the ODT in the experiments: the shifting up of the ODT and direct transitions between disorder and the ordered phases.

Scientific Motivation

1. Use Monte Carlo simulations (more details in a previous project) to investigate the effects of fluctuations on diblock copolymer melts.
2. Provide a real-time, visual output for the underlying Monte Carlo model to facilitate quick inspection of the system configuration/properties and aid parameter tuning.

Technical Motivation

1. The fast inspection of system configurations and properties aids in quickly tuning parameters that can then be used in faster simulations (no front-end) for longer runs.
2. Thermodynamic averages/statistics can be collected and displayed in charts in real-time as the simulation progresses (quicker than having to plot data files from a command-line program).

Broad Simulation Requirements

1. Create a visual front-end that allows the user to specify input parameters and displays the system configuration/data in real-time.
2. Define a 3D lattice and generate an intial configuration.
3. Use a Monte Carlo model (treats fluctuation effects) as the underlying engine on the back-end.
4. Display performance-related information, such as time taken per Monte Carlo update and simulation progress.
5. Display quantities such as the structure function, $S(k)$ and energy, $U$.
6. Allow the user to control how often the configuration display is updated (frequent updates will limit simulation speed).

Technical Simulation Requirements

1. Employ the lattice Monte Carlo model detailed previously (but without parallel tempering) for the back end.
2. Use C# language as it is easy to build an interactive front-end with Windows Presentation Foundation (WPF).
3. Port the Monte Carlo code from C/C++ to C# (relatively simple due to the similarity of the languages).
4. Run the Monte Carlo simulation on background worker threads to keep the front-end responsive.
5. Use the Helix Toolkit library (allows rotation and zooming) to display the system configuration.

In statistical physics, a system with fixed volume and number of particles in contact with a heat bath at a fixed temperature is referred to as a canonical ensemble. The partition function for a canonical ensemble, such as our diblock copolymer melt, is given by $$Z=\sum _i{e}^{-\beta E_i},$$ where $E_i$ is the energy of microstate (specific configuration) $i$ and $\beta =1/k_{B}T$. Using the partition function, it is possible to calculate the thermodynamic average of a system observable, $A$, via: $$⟨A⟩=\frac{1}{Z}\sum _i{A}_i{e}^{-\beta E_i},$$ where $A_i$ is the value of observable $A$ in state $i$. For example, the thermodynamic average of the energy would be: $$⟨E⟩=\frac{1}{Z}\sum _i{E}_i{e}^{-\beta E_i}.$$ The problem is that the number of states in a system becomes so large that we can't possibly count them and calculate the partition function directly. This is where Monte Carlo simulations come in...

Monte Carlo simulations create what is known as a 'Markov chain', where the next state to be visited depends only on the current state. The parition function of the system is effectively sampled by visiting each system-state with the correct probability - less-likely states are visited fewer times that likely states. Samples of observables (such as energy) are taken periodically as the system evolves and subsequently averaged. The longer the simulations runs for (more states visited), the more accurately the Monte Carlo averages will approximate the true thermodynamic values. The algorithm used to update the system in this project is summarised as follows:

Monte Carlo Metropolis Algorithm

1. The system begins in a state, $i$, with energy $U_i$.
2. By performing one of the prescribed Monte Carlo moves, a small perturbation is made to the initial configuration in an attempt to move the system to state $j$, with energy $U_j$.
3. The difference in energy, $\mathrm{\Delta }U=U_j-Ui$, that would result from moving the system from state $i\to j$ is calculated.
4. If $\mathrm{\Delta }U\le 0$, then the accepted move from state $i\to j$ is accepted, the system configuration is updated and we move to step 6.
5. If $\mathrm{\Delta }U>0$, the system is attempting to adopt a higher energy state. In this case, generate a random number, $R$. If $R<\exp (-\mathrm{\Delta }U/k_{B}T)$ then the move from $i\to j$ is accepted and the system configuration updated. Otherwise, the system remains in the initial state, $i$.
6. Take samples contributing to Monte Carlo averages.
7. Go to step 1 and repeat.

Benefits of Monte Carlo

There are many benefits to using the Monte Carlo method in this project, but the particularly significant ones are as follows:

Lattice Model Description

Each molecule in the system has $N_A$ A-type monomers and $N_B$ B-type molecules, such that there are $N=N_A+N_B$ in total in a particular diblock. The monomers are placed on a lattice such that each lattice site may be vacant or occupied by a single monomer and bonded monomers must occupy nearest-neighbour sites. Although it would be more physically realistic to perform the simulations in free-space, a lattice construction is far less computationally expensive and allows for more Monte Carlo statistics to be collected per unit time.

The lattice is contructed by starting with an $L\times L\times L$ cubic lattice with $V_c=L^3$ sites and lattice constant, $d$. The allowed lattice vectors can be described by ${\mathbf{r}}_i=d(h,k,l)$, where $i=1,2,\dots ,L^3$ is the site number and $h,k,l=0,1,\dots ,L-1$. By allowing only sites where $h+k+l$ is even, a face-centred cubic (fcc) lattice with bond length $b=\sqrt{2}d$ and $V=V_c/2=L^3/2$ sites results. An fcc construction is used as its large coordination number (number of nearest-neighbours a particular site has) of $z=12$ helps to alleviate the effects of the lattice. Periodic boundary conditions are implemented in all three dimensions to reduce the effects of the finite lattice size, requiring that $L$ is an even number. The periodic boundary conditions mean that a molecule can travel off a face of the lattice and continue on the opposite face.

To allow the molecules room to move around and for the Monte Carlo algorithm to become efficient, the lattice is filled to a copolymer concentration of ${\varphi }_c\equiv nN/V\approx 0.8$, where n is the total number of molecules in the system. Ideally, all lattice sites would be occupied, but this would make the Monte Carlo algorithm more complex, reducing its speed and the size of the systems that could be investigated.

Molecular interactions within the system occur only between nearest-neighbour A and B monomers. The strength of the interaction is given by ${ϵ}_{AB}$, from which a dimensionless Flory-Huggins parameter $$\chi \equiv \frac{z{ϵ}_{AB}}{k_{B}T}$$ is defined.

At each Monte Carlo step (MCS), an attempt is made to change the system slightly by performing one of the available local moves. If the attempted move results in double occupancy of any lattice site, it is immediately rejected according to the Metropolis algorithm described in section 'Monte Carlo Methods'.

Local Monte Carlo Moves

Selection and implementation of the allowed Monte Carlo moves can have a large effect of the efficiency of the Monte Carlo algorithm. The moves should quickly update the system's configuration and efficiently sample its phase space. At each MCS, one of three moves: snake, crank or flip (figs 7a-c), is randomly attempted according to the predefined probabilities: $P_{snake}=0.45$, $P_{crank}=0.45$ and $P_{flip}=0.10$.

Setup and Input Parameters

After running the program, the first window presented to the user is shown in Figure 6. The "Input Parameters" section is populated with some initial values that can be edited by the user (the individual parameters are described in the next section).

Once finished editing the parameters, clicking the "Set up" button will validate that the parameters represent a valid simulation. If validation fails, the user is represented with an error message instructing them on how to fix the issue(s). After passing validation, the Monte Carlo back-end proceeds to generate an initial configuration as shown in Figure 7.

Setup and Input Parameters

$L$
Number of sites along the $x,y,z$ dimensions of the simulation box.

$M$
Total number of polymers in the system.

$N$
Total number of monomers in an $AB$-diblock copolymer.

$N_A$
Number of monomers in the $A$-block of an $AB$-diblock copolymer.

$\chi N$
The Flory-Huggins interaction parameter multipled by $N$. Dictates the energy penalty of $A-B$ monomer interactions.

$n_{\text{relax}}$
Length of relaxation period before sampling system statistics (in units of Monte Carlo steps per monomer).

$n_{\text{estim}}$
Length of statistics gathering period (in units of Monte Carlo steps per monomer).

$f_{\text{sample}}$
Period between taking samples in the estimation phase (in units of Monte Carlo steps per monomer).

$(P_{\text{snake}},P_{\text{flip}},P_{\text{crank}})$
Cumulative probabilities with which the snake, flip and crank moves are attempted (must sum to one). E.g., $(0.45,0.55,1.0)$ attempts the snake move 45% of the time, the flip move 10% of the time and the crank move the remaining 45% of the time.

Initial Configuration and Interacting with the Viewer

Once the input parameters have been validated, the Monte Carlo model sets up an initial configuration that is subsequntly displayed in the configuration viewer by the front-end (Figure 7). The user is then able to interact with the configuration:

Rotation:
Right-click and hold on the configuration. Moving the mouse then allows the configuration to be rotated in any direction.

Zoom:
Scroll the mouse wheel on the configuration to zoom in/out.

Translate:
Push and hold the middle mouse button on the configuration. Moving the mouse then allows the configuration to be translated within the viewer.

Reset View:
Double clicking the middle mouse button in the configuration viewer will reset the view to its default position.

Running the Simulation and Viewing its Status

After inspecting the initial (non-equilibrium) configuration, the user can click the "Run" button to start the Monte Carlo simulation. The simulation's progress is displayed to the user in the "Status" Section, as shown in Figure 8 and decribed below:

Equilibration:
The percentage of the simulation's equilibration (relaxation) period completed. Displayed as text and a progress bar.

Statistics:
The percentage of the simulation's statistics period completed. Displayed as text and a progress bar.

Sweep #:
The step number reached by the simulation at the last update of the front-end (in units of Monte Carlo steps per monomer).

ms/sweep:
The average amount of time (in milliseconds) for the system to attempt an update of every monomer in the system.

Using the Controls

While the simulation is running, the "Controls" section allows the user to alter how the configuration is displayed, as well as how often. The various controls are detailed as followed:

"Update Plot" Checkbox:
Switches on/off updates of the configuration display (the simulation continues in the background).

"AB-A-B" Radio Group:
Selects whether to display all monomers (AB), A-monomers (A) or B-monomers (B).

"Update Freq" Slider:
Specifies how often to update the configuration displayed on the front-end.

"Shell Depth" Slider:
The number of lattice layers (from the outside in) that are displayed when drawing the configuration.

Efficiency Considerations

It should be noted that displaying the configuration too often can negatively impact the speed of the Monte Carlo simulation. The size of the effect depends upon a number of factors, including simulation box size, but can be mitigated to a large extent by our choice of display control parameters:

1. We are usually interested in observing the faces of the simulation box, which obscure the inside from view. Using the "Shell Depth" slider to display only one or two layers stops the internal (unobservable) layers from being drawn, saving computational effort and increasing simulation speed. If the internal layers need to be observed (e.g., checking domain purity while displaying only A or B monomers), the shell depth can be momentarily increased.
   
   
2. Reducing how often the configuration is redrawn via the "Update Freq" slider will take fewer computational resources away from the Monte Carlo simulation working in the background. Alternatively, unticking the "Update Plot" checkbox will have the most drastic effect.

When altering the display parameters, the user should observe their effects on the "ms/sweep" in the "Status" section (a lower number is better).

Coming soon -> This project is still a work in progress.

Although the core simulation and front-end are functioning correctly, there are a number of features that are not yet fully implemented, including real-time chart updates that will display the system's thermodynamic averages (concept already proven). As such, the program will be made available once these features have been added and tested.