Variable Cell Formalism

This derivation of this term can be found in Barrat2005, Villet2014 Appendix B, in unpublished notes by Kris Delaney or in Josh’s notes NB#2 p140-148. Here we follow the notation of Kris Delaney.

Note: add bold math using pmb{}

Following Barrat2005 and GHF section 5.3.5 (p 252) we can conduct simulations in the constant tension $N\tau P$ ensemble. This ensemble is particularly useful for calculating the cell size of periodic mesophases, which can be found by adjusting the cell to minimize the internal stresses.

For a cell with box vectors $e{e}_1$, $e{e}_2$, $e{e}_3$, we define the cell shape tensor

$$hh=\left(\begin{array}{ccc}e{e}_1& e{e}_2& e{e}_3\end{array}\right)$$

The metric tensor is defined as

$$gg\equiv h{h}^{T}hh$$

To find the optimum box size, we seek to minimize $\frac{\partial H_V}{\partial hh}$ where $H_V\equiv H/V$ is the intensive energy density. This derivative is related to the internal stress operator $V\sigma {\sigma }_{int}=\frac{\partial H}{\partial hh}h{h}^T$. For a general multi-chain Hamiltonian, $H\sim H_{loc}-\sum _p\frac{C\tilde{V}{\varphi }_p}{{\alpha }_p}\ln Q_p$. The stress operator is thus

$$V\sigma {\sigma }_{int}=\frac{\partial H_{loc}}{\partial hh}h{h}^T-\sum _p\frac{CV{\varphi }_p}{{\alpha }_p{Q}_p}\frac{\partial Q_p}{\partial hh}h{h}^T$$

For most of the polymer models described in GHF (i.e. without Coulomb interactions), $\frac{\partial H_{loc}}{\partial hh}=0$ and the calculation of $V\sigma {\sigma }_{int}$ involves the calculation of $\frac{\partial Q_p}{\partial hh}$. To derive this expression, one first shows that $\frac{\partial Q_p}{\partial hh}=2hh\frac{\partial Q_p}{\partial gg}$. Recalling that $Q_p=\int dxx{q}_p(xx,{\alpha }_p)=\frac{1}{V}\int drr{q}_p(rr,{\alpha }_p)$ where $q_p$ satisfies the general PDE $\frac{\partial }{\partial s}q(xx,s)=\mathcal{L}q(xx,s)$, one can derive

$$\frac{\partial Q_p}{\partial gg}=\int dxx{\int }_0^{{\alpha }_p}{q}_p^{†}(xx,{\alpha }_p-s)\frac{\partial \mathcal{L}}{\partial gg}{q}_p(xx,s).$$

For a continuous Gaussian chain model, $\mathcal{L}=\mathrm{\nabla }\mathrm{\nabla }\mathrm{\nabla }\mathrm{\nabla }-iw(xx)=g{g}^{-1}:\mathrm{\nabla }{\mathrm{\nabla }}_x\mathrm{\nabla }{\mathrm{\nabla }}_x-iw(xx)$ and so

$$\frac{\partial Q_p}{\partial gg}=-\int dxx{\int }_0^{{\alpha }_p}{q}_p^{†}(xx,{\alpha }_p-s)h{h}^{-1}\mathrm{\nabla }\mathrm{\nabla }\mathrm{\nabla }{\mathrm{\nabla }}^{T}h{h}^{-T}{q}_p(xx,s).$$

Combining all terms, the final result for an unsmeared model is

$$V\sigma {\sigma }_{int}=\frac{\partial H_{loc}}{\partial hh}h{h}^T+\sum _p\frac{2CV{\varphi }_p}{{\alpha }_p{Q}_p}\frac{1}{V}\int drr{\int }_0^{{\alpha }_p}ds{q}_p^{†}(rr,{\alpha }_p-s)\mathrm{\nabla }\mathrm{\nabla }\mathrm{\nabla }{\mathrm{\nabla }}^T{q}_p(rr,s)$$

In practice, the arguement of the $\int ds$ integral can be efficiently computed in k-space

$$V\sigma {\sigma }_{int}=\sum _p\frac{2CV{\varphi }_p}{{\alpha }_p{Q}_p}\int drr{\int }_0^{{\alpha }_p}ds\sum _{k}k{q}_p^{†}(-kk,{\alpha }_p-s)(-kkk{k}^T)q_p(kk,s)$$

where we have now set $\frac{\partial H_{loc}}{\partial hh}=0$.

For smeared models where $Q_p[\mathrm{\Gamma }\ast w]$, an addition term emerges since $\mathrm{\Gamma }$ depends on $gg$. For smeared models the stress is

$$\begin{array}{rl}V\sigma {\sigma }_{int}=& \sum _p\frac{2CV{\varphi }_p}{{\alpha }_p{Q}_p}\int drr{\int }_0^{{\alpha }_p}ds\sum _{k}k{q}_p^{†}(-kk,{\alpha }_p-s)(-kkk{k}^T)q_p(kk,s)\\ & +\sum _{p}CV\sum _{i=1}^S{a}_i\sum _{kk}{\hat{\rho }}_i(-kk)\hat{\mathrm{\Gamma }}(kk){\hat{w}}_i(kk)(kkk{k}^T)\end{array}$$

where the sum over i runs over each species, $a_i$ is the smearing range, and ${\rho }_i$ is the the density operator,

$${\rho }_i(xx)=\frac{{\varphi }_p}{{\alpha }_p{Q}_p}{\int }_{s\in i}ds{q}^{†}(xx,{\alpha }_p-s)q(xx,s).$$

Once ${\sigma }_{int}$ is obtained, $H_V$ can be minimized by updating $hh$ using a forward Euler scheme

$$\begin{array}{rl}h{h}^{j+1}& =h{h}^j-\delta t{\lambda }_h{\left[\frac{\partial H_V}{\partial hh}\right]}^j\end{array}$$

where $\delta t$ is the timestep used to update the $\mu$ fields, ${\lambda }_h$ adjusts the rate of the cell updates relative to the fields, and $\frac{\partial H_V}{\partial hh}={\sigma }_{int}h{h}^{-T}$ (as defined above).

LEGACY DOCUMENTATION

(This corresponds to notation from GHF and Barrat2005, above I use notation from Kris Delaney)

The only term in $H_V$ that depends on $hh$ is the $\log Q$ term. For these models

$$\frac{\partial H_V}{\partial hh}=\beta hh\mathrm{\Sigma }\mathrm{\Sigma }$$

where $\mathrm{\Sigma }$ is a symmetric tensor that for molecule p whose elements ${\mathrm{\Sigma }}_{p,\alpha \beta }$ are given by

$$\beta {\mathrm{\Sigma }}_{p,\alpha \beta }=\frac{n_p}{3VQ}\int dxx{\int }_0^{N_p}ds[b(s){]}^2{g}_{\gamma \alpha }^{-1}q(xx,s)\frac{{\partial }^2}{\partial x_{\gamma }\partial x_{\delta }}{q}^{†}(xx,s)g_{\delta \beta }^{-1}$$

where the sum over repeated indicies is implied (i.e. Einstein notation).

Note

The prefactor in Barrat2005 Eq 36 is $\frac{n}{V}\frac{2R_{g,0}^2}{Q}$ and the $ds$ integral runs from 0 to 1. This comes from how the units are non-dimentionalized, but I never got to the bottom of what the precise prefactor should be. There’s a chance that the precise magnitude of $\mathrm{\Sigma }$ is incorrect in the code. However, since we’re usually interested in where $\mathrm{\Sigma }=0$, it shouldn’t matter (at least for the moment).

Since ${\mathrm{\Sigma }}_p$ involves $q(xx,s)$ and $q^{†}(xx,s)$ it is straightforward to calculate ${\mathrm{\Sigma }}_p$ while solving the modified diffusion equation for the propagators.

Following a computation trick described by Kris Delaney, $\mathrm{\Sigma }$ can efficiently be substituting the definition of the fourier transform of $q(rr,s)=\frac{1}{V}\sum _{k}k\hat{q}(kk,s)e^{ikk\cdot xx}$. Also recalling that $\frac{{\partial }^2}{\partial x_{\gamma }\partial x_{\delta }}\to -k_{\gamma }{k}_{\delta }$ in fourier space,

$$\beta {\mathrm{\Sigma }}_{p,\alpha \beta }=\frac{n_p}{3V{Q}_p}{\int }_0^{N_p}ds[b(s){]}^2\frac{1}{V^2}\sum _{k}k\hat{q}(kk,s)g_{\gamma \alpha }^{-1}(-k_{\gamma }{k}_{\delta })g_{\delta \beta }^{-1}{\hat{q}}^{†}(kk,s)$$

where we have used the identity $\int dxx{e}^{+i(kk+k{k}^{\prime })\cdot xx}=\delta (kk+k{k}^{\prime })$. The term inside the $ds$ integral is computed for each s during the solution of the propagator. Also note that $g{g}^{-1}(-kkkk)g{g}^{-1}$ only needs to be computed once for each call shape $hh$.

By switching to the Duchs2014 notation described in model_chi_ab.rst for a system of many polymers

$$\beta \mathrm{\Sigma }\mathrm{\Sigma }=\sum _p\frac{C{\varphi }_p}{{\alpha }_p}\mathrm{\Sigma }{\mathrm{\Sigma }}_p$$

where in this equation the first $\sum _p$ (i.e. \sum_p) is the sum over each polymer p and the second $\mathrm{\Sigma }{\mathrm{\Sigma }}_p$ (i.e. \Sigma_p) is a symmetric tensor related to the stress (yes, I realize this really isn’t the best notation).

Once $\mathrm{\Sigma }$ is obtained, $H_V$ can be minimized by updating $hh$ using a forward Euler scheme

$$h{h}^{j+1}=h{h}^j-\delta t{\lambda }_h{\left[\frac{\partial H_V}{\partial hh}\right]}^j$$

where $\delta t$ is the timestep used to update the $\mu$ fields, ${\lambda }_h$ adjusts the rate of the cell updates relative to the fields, and $\frac{\partial H_V}{\partial hh}=\beta hh\mathrm{\Sigma }\mathrm{\Sigma }$ (as defined above).

Finally, we define the “internal stress operator”

$$\tilde{\sigma }\equiv hh\mathrm{\Sigma }h{h}^T$$