Problem: given a series molecules $i=1,\dots,N$, connected with junction atoms $\mathbf{x}_i^k$, for $k=1,2,…m$ if there are $m_i$ atoms from molecule $i$ that are bonded with other molecules. The connection bonds are defined as $\mathbf{b}_i := \lbrace\mathbf{x}_j^{k_1}-\mathbf{x}_i^{k_2}|\forall k_1, k_2\, \mathrm{and}\, \forall j\in \mathrm{neighbor\, of}\, i\rbrace$, assume that the center-of-mass of the molecules are fixed, but are rotatable, minimize $\sum_i |\mathbf{b}_i|^2$ with a series of rotation matrices $R_i$:

$$ R_i = \min_{R_i}(L:= \sum_i \sum_{j \in \mathrm{neighbor\, of}\, i} \sum_{k_1,k_2} |\mathrm{pbc}(R_j \mathbf{x}_{j}^{k_1} – R_i \mathbf{x}_{i}^{k_2})|^2) $$

with constraint

$$ c=\sum_i |R_i^TR_i – I|^2=0 $$

the constraint ensures $R^T R =I$ for each rotation matrix thus $R_i$ are rotation matrices. The periodic boundary condition is introduced for most simulation systems.

Minimization procedure includes calculation of jacobian of loss and constraint functions:

$$ \frac{\partial L}{\partial R_i} = -2 \sum_{j \in \mathrm{neighbor\, of}\, i} \sum_{k_1, k_2} \mathrm{pbc}(R_{j} \mathbf{x}_{j}^{k_1} – R_i \mathbf{x}_i^{k_2})^T \mathbf{x}_i^{k_2} $$

and

$$\frac{\partial c}{\partial R_i} = 4 (R_i R_i^T – I) R_i$$

The structure factor, commonly denoted as $S(\mathbf{q})$, is an essential characterization quantity in molecular simulations. Its calculation is straightforward based on its definition. Since periodic boundary conditions are typically used in these simulations, we have

$$S(\mathbf{q}) := \mathcal{FT}\lbrace\langle\rho(\mathbf{r})\rho(\mathbf{0})\rangle\rbrace(\mathbf{q})$$

which actually represents a circular convolution. This means that using the Fast Fourier Transform (FFT) can dramatically accelerate the computation compared to performing the calculation directly.

The FFT-based procedure proceeds as follows:

1. Compute the density field $\rho(\mathbf{r})$, which is a sum of Dirac delta functions approximated by a histogram.
2. Use the FFT to calculate the Fourier transform $\hat{\rho}(\mathbf{q})$.
3. Obtain the structure factor by taking the modulus squared of the transform:
   $$ S(\mathbf{q}) = |\hat{\rho}(\mathbf{q})|^2. $$
4. Average over the modulus to yield
   $$ S(q) = \frac{\int S(\mathbf{q})\delta(|\mathbf{q}|-q)\,\mathrm{d}\mathbf{q}}{\int\delta(|\mathbf{q}|-q)\,\mathrm{d}\mathbf{q}}. $$

In terms of efficiency, if the simulation box is divided into $N$ bins (for example, in 3D, $N = N_x N_y N_z$), the FFT has a computational complexity of $O(N\log(N))$, and the averaging step is $O(N)$. By contrast, a direct method would involve a loop over all particles and $N$ bins in $\mathbf{q}$ space—and typically, the number of particles is much larger than $\log(N)$. Moreover, because most simulation data are real-valued, employing an FFT designed for real inputs (rFFT) can further improve performance.

Below is an example code for a 3D case, where the indices i, j, and k correspond to $N_x$, $N_y$, and $N_z$, respectively.

fftn(a) == np.concatenate([rfftn(a), conj(rfftn(a))[-np.arange(i),-np.arange(j),np.arange(k-k//2-1,0,-1)]], axis=-1)

Accuracy Considerations: Effect of Bin Size

The histogram bin size should be smaller than half the minimum pairwise distance in the system, in accordance with the Nyquist-Shannon sampling theorem. For instance, in Lennard-Jones systems, a bin size of approximately $0.4\sigma$ is advisable. However, if you are only interested in small $q$ values—as is often the case in self-assembly systems—full pairwise resolution is not necessary; a sampling interval that is just smaller than half of the domain of interest may suffice.

Randomness Within Bins

Particles are assumed to be randomly distributed within each bin, especially when employing large bins. The Fourier transform of the particle positions is given by

$$\sum_i \exp(-\mathrm{i}\mathbf{q}\cdot\mathbf{r}_i)$$

This sum can be decomposed as

$$\sum_\mathrm{bin} \exp(-\mathrm{i}\mathbf{q}\cdot\mathbf{r}_\mathrm{bin})\left(\sum_i \exp(-\mathrm{i}\mathbf{q}\cdot\delta\mathbf{r}_i)\right)$$

where each particle’s position is written as

$$\mathbf{r}_i = \mathbf{r}_\mathrm{bin} + \delta\mathbf{r}_i$$

Assuming the displacements $\delta\mathbf{r}_i$ are randomly distributed within the bin, the inner sum can be approximated by

$$n_\mathrm{bin}\langle\exp(-\mathrm{i}\mathbf{q}\cdot\delta\mathbf{r})\rangle$$

with the average taken over the distribution of $\delta\mathbf{r}$. For a one-dimensional example, if $\delta r$ is uniformly distributed in the interval $(-L/2, L/2)$ (where $L$ is the bin size), the average is given by $\mathrm{sinc}(qL/2)$. Multiplying by this sinc correction during step 4 yields a properly corrected $S(q)$.

This strategy ensures that both efficiency and accuracy are maintained when computing the structure factor in molecular simulations.

In my previous post, I developed a GPC post-processing program that uses the continuous wavelet transform (CWT) to detect and correct baselines, and the MH equation to convert GPC data into a molecular weight distribution. In this follow-up work, I determine the molecular weight distributions in a two-step synthesis process. In the first step, polymers are pre-polymerized with a distribution denoted as $P_1(M)$. In the second step, these polymers continue to grow, ultimately yielding a final chain-length distribution $Q(M)$. The goal here is to extract the molecular weight distribution for the second step, denoted as $P_2(M)$.

I present a straightforward method to calculate the molecular weight distribution of the second block using the GPC data from the preceding block and the final diblock copolymer. Assuming that we already know the distribution of the preceding block, $P_1(n)$, it is generally expected that the growth of the second block depends on the preceding chain length $n$. Consequently, the joint probability density function (pdf) for the two blocks is given by
$$ Q(n, m) = P_1(n)P_2(m, n)$$
and the pdf of the final diblock copolymer is obtained as
$$ Q(x) = \int Q(n, x-n) \mathrm{d}n $$

It is reasonable to assume that $P_2(m, n)$ can be factorized into the product of two functions, say $f(m)g(n)$, by neglecting higher-order correlations between $m$ and $n$. In this factorization, $g(n)$ acts as a scaling factor on the distribution $P_1(n)$, and the calculation of $P(x)$ becomes a convolution. The simplest case is when $g(n)=1$, meaning that the growth of the second block is independent of the length of the preceding block. Alternatively, one might assume $g(n) \sim n^{-1}$ to account for diffusion effects, implying that shorter preceding chains tend to grow a longer second block.

The evaluation of the pdf for the second block proceeds in three steps:

1. Determine the range of molecular weights as $(x_{min} – n_{max},\, x_{max} – n_{min})$, where $x$ is the chain length of the final diblock copolymer and $n$ is the length of the preceding block.
2. Interpolate the GPC-derived distributions onto an evenly spaced molecular weight grid, setting any negative values to zero.
3. Deconvolute $Q(x)$ with $P_1(n)$.

This method provides a straightforward pathway to isolate the molecular weight distribution of the second synthetic block from the overall diblock copolymer distribution.

Consider a polymer chain dissolved in a solvent, where its behavior is influenced by the Flory parameter, denoted by $\chi$. As the interaction parameter $\chi$ varies, the chain undergoes several conformational transitions:

• For $\chi > \chi_c$, the polymer collapses into a dense, globular structure with a radius of gyration scaling as $R_g \sim N^{1/3}$.

• At the critical point—commonly known as the $\theta$ point—where $\chi \sim \chi_c$, the chain behaves like an ideal random coil, characterized by $R_g \sim N^{1/2}$.

• In a good solvent, when $\chi < \chi_c$, the polymer chain is swollen and stretched, with $R_g \sim N^{3/5}$.

A “universal” expression that captures how the radius of gyration $R_g$ changes with $\chi$ is given by

  $R_g(\chi) = R_g^\mathrm{glob} + \frac{R_g^\mathrm{coil}-R_g^\mathrm{glob}}{1 + \exp[(\chi-\chi_c)/\Delta \chi]}$,

where $\Delta \chi$ quantifies the width of the conformational transition region and $\chi_c$ locates the transition.

To estimate the window width $\Delta \chi$, we first define an order parameter that measures deviations of the polymer’s size from its value at the $\theta$ point. We write

  $R = R_0 (1 + m)$,

with $R_0 \sim N^{1/2}b$, where $b$ is the monomer size and $m$ represents the fractional deviation from the ideal size.

Next, we expand the free energy $F(m)$ in powers of $m$. Close to the $\theta$ point, the expansion takes the form

  $\frac{F(m)}{k_BT} \simeq C_1 N (\chi-\chi_c) m^2 + C_2 N m^4 + \cdots$,

where $C_1$ and $C_2$ are constants, and $k_BT$ is the thermal energy.

In the mean-field (infinite-system) limit, the phase transition is sharp. However, for finite systems, fluctuations smear the transition, leading to a rounded behavior characterized by a finite width $\Delta \chi$ defined by $|\chi-\chi_c|$.

At the edge of the transition, the quadratic and quartic terms in the free energy become comparable. Equating these contributions gives

  $C_1 N\,\Delta\chi\,m^2 \sim C_2 N\,m^4,$

which implies

  $m^2 \sim \frac{C_1}{C_2}\,\Delta\chi$.

Furthermore, the transition is significantly broadened when the free energy barrier, estimated by the quartic term, is of the order of unity (in units of $k_BT$). That is, when

  $N C_2\, m^4 \sim O(1).$

Substituting our earlier estimate for $m^2$ into this condition, we obtain

  $N C_2 \left(\frac{C_1}{C_2}\Delta\chi\right)^2 \sim O(1),$

which simplifies to

  $\Delta\chi^2 \sim \frac{1}{N}.$

Thus, the rounding of the transition in terms of $\chi$ scales as

  $\Delta\chi \sim \frac{1}{\sqrt{N}}.$

Technical Approach:

Preparation of Source Images

• Capture Three Photos:
  • Full Moon: Take a clear photograph of the full moon.
  • Two Waning Moons: Photograph the waning moon twice.
    • Overexposed Waning Moon: For one of the waning moon images, apply an exposure compensation of +3 levels to create an overexposed effect, which will enhance the halo effect around the moon.

(more…)