Variable Cell =================== 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 :math:`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 :math:`\pmb{e}_1`, :math:`\pmb{e}_2`, :math:`\pmb{e}_3`, we define the cell shape tensor .. math:: \pmb{h} = \begin{pmatrix} \pmb{e}_1 & \pmb{e}_2 & \pmb{e}_3 \end{pmatrix} The metric tensor is defined as .. math:: \pmb{g} \equiv \pmb{h}^T \pmb{h} To find the optimum box size, we seek to minimize :math:`\frac{\partial H_V}{\partial \pmb{h}}` where :math:`H_V \equiv H / V` is the *intensive* energy density. This derivative is related to the internal stress operator :math:`V\pmb{\sigma}_{int} = \frac{\partial H}{\partial \pmb{h}} \pmb{h}^T`. For a general multi-chain Hamiltonian, :math:`H \sim H_{loc} - \sum_p \frac{ C \tilde{V} \phi_p}{\alpha_p} \ln Q_p`. The stress operator is thus .. math:: V\pmb{\sigma}_{int} = \frac{\partial H_{loc}}{\partial \pmb{h}} \pmb {h}^T - \sum_p \frac{ C V \phi_p}{\alpha_p Q_p} \frac{\partial Q_p}{\partial \pmb{h}} \pmb{h}^T For most of the polymer models described in GHF (i.e. without Coulomb interactions), :math:`\frac{\partial H_{loc}}{\partial \pmb{h}} = 0` and the calculation of :math:`V\pmb{\sigma}_{int}` involves the calculation of :math:`\frac{\partial Q_p}{\partial \pmb{h}}`. To derive this expression, one first shows that :math:`\frac{\partial Q_p}{\partial \pmb{h}} = 2 \pmb{h} \frac{\partial Q_p}{\partial \pmb{g}}`. Recalling that :math:`Q_p=\int d\pmb{x} q_p(\pmb{x},\alpha_p)=\frac{1}{V} \int d\pmb{r} q_p(\pmb{r},\alpha_p)` where :math:`q_p` satisfies the general PDE :math:`\frac{\partial}{\partial s} q(\pmb{x},s) = \mathcal{L} q(\pmb{x},s)`, one can derive .. math:: \frac{\partial Q_p}{\partial \pmb{g}} = \int d\pmb{x} \int_0^{\alpha_p} q_p^{\dagger}(\pmb{x},\alpha_p - s) \frac{\partial \mathcal{L}}{\partial \pmb{g}} q_p(\pmb{x},s). For a continuous Gaussian chain model, :math:`\mathcal{L} = \pmb{\nabla}\pmb{\nabla} - i w(\pmb{x}) = \pmb{g}^{-1} : \pmb{\nabla}_x \pmb{\nabla}_x - i w(\pmb{x})` and so .. math:: \frac{\partial Q_p}{\partial \pmb{g}} = - \int d\pmb{x} \int_0^{\alpha_p} q_p^{\dagger}(\pmb{x},\alpha_p - s) \pmb{h}^{-1} \pmb{\nabla} \pmb{\nabla}^T \pmb{h}^{-T} q_p(\pmb{x},s). Combining all terms, the final result for an unsmeared model is .. math:: V\pmb{\sigma}_{int} = \frac{\partial H_{loc}}{\partial \pmb{h}} \pmb {h}^T + \sum_p \frac{ 2 C V \phi_p}{\alpha_p Q_p} \frac{1}{V}\int d\pmb{r} \int_0^{\alpha_p} ds q_p^{\dagger}(\pmb{r},\alpha_p-s) \pmb{\nabla}\pmb{\nabla}^T q_p(\pmb{r},s) In practice, the arguement of the :math:`\int ds` integral can be efficiently computed in k-space .. math:: V\pmb{\sigma}_{int} = \sum_p \frac{ 2 C V \phi_p}{\alpha_p Q_p} \int d\pmb{r} \int_0^{\alpha_p} ds \sum_\pmb{k} q_p^{\dagger}(-\pmb{k},\alpha_p-s) (-\pmb{k}\pmb{k}^T) q_p(\pmb{k},s) where we have now set :math:`\frac{\partial H_{loc}}{\partial \pmb{h}} = 0`. For smeared models where :math:`Q_p[\Gamma * w]`, an addition term emerges since :math:`\Gamma` depends on :math:`\pmb{g}`. For smeared models the stress is .. math:: V\pmb{\sigma}_{int} = &\sum_p \frac{ 2 C V \phi_p}{\alpha_p Q_p} \int d\pmb{r} \int_0^{\alpha_p} ds \sum_\pmb{k} q_p^{\dagger}(-\pmb{k},\alpha_p-s) (-\pmb{k}\pmb{k}^T) q_p(\pmb{k},s) \\ &+ \sum_p C V \sum_{i=1}^S a_i \sum_{\pmb{k}} \hat{\rho}_i(-\pmb{k}) \hat{\Gamma}(\pmb{k}) \hat{w}_i(\pmb{k}) (\pmb{k}\pmb{k}^T) where the sum over `i` runs over each species, :math:`a_i` is the smearing range, and :math:`\rho_i` is the the density operator, .. math:: \rho_i(\pmb{x}) = \frac{\phi_p}{\alpha_p Q_p} \int_{s \in i} ds q^{\dagger}(\pmb{x}, \alpha_p - s) q(\pmb{x},s). Once :math:`\sigma_{int}` is obtained, :math:`H_V` can be minimized by updating :math:`\pmb{h}` using a forward Euler scheme .. math:: \pmb{h}^{j+1} &= \pmb{h}^j - \delta t \lambda_h \left[ \frac{\partial H_V}{\partial \pmb{h}} \right] ^j \\ where :math:`\delta t` is the timestep used to update the :math:`\mu` fields, :math:`\lambda_h` adjusts the rate of the cell updates relative to the fields, and :math:`\frac{\partial H_V}{\partial \pmb{h}} = \sigma_{int} \pmb{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 :math:`H_V` that depends on :math:`\pmb{h}` is the :math:`\log Q` term. For these models .. math:: \frac{\partial H_V}{\partial \pmb{h}} = \beta \pmb{h} \pmb{\Sigma} where :math:`\Sigma` is a symmetric tensor that for molecule `p` whose elements :math:`\Sigma_{p,\alpha\beta}` are given by .. math:: \beta \Sigma_{p,\alpha\beta} = \frac{n_p}{3VQ} \int d\pmb{x} \int_0^{N_p} ds [b(s)]^2 g^{-1}_{\gamma\alpha} q(\pmb{x},s) \frac{\partial^2}{\partial x_{\gamma} \partial x_{\delta}} q^\dagger(\pmb{x},s) g^{-1}_{\delta\beta} where the sum over repeated indicies is implied (i.e. Einstein notation). .. note:: The prefactor in Barrat2005 Eq 36 is :math:`\frac{n}{V}\frac{2 R_{g,0}^2}{Q}` and the :math:`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 :math:`\Sigma` is incorrect in the code. However, since we're usually interested in where :math:`\Sigma=0`, it shouldn't matter (at least for the moment). Since :math:`\Sigma_p` involves :math:`q(\pmb{x},s)` and :math:`q^\dagger(\pmb{x},s)` it is straightforward to calculate :math:`\Sigma_p` while solving the modified diffusion equation for the propagators. Following a computation trick described by Kris Delaney, :math:`\Sigma` can efficiently be substituting the definition of the fourier transform of :math:`q(\pmb{r},s) = \frac{1}{V} \sum_\pmb{k} \hat{q}(\pmb{k},s) e^{i\pmb{k}\cdot\pmb{x}}`. Also recalling that :math:`\frac{\partial^2}{\partial x_\gamma \partial x_\delta} \rightarrow -k_\gamma k_\delta` in fourier space, .. math:: \beta \Sigma_{p,\alpha\beta} = \frac{n_p}{3VQ_p} \int_0^{N_p} ds [b(s)]^2 \frac{1}{V^2} \sum_\pmb{k} \hat{q}(\pmb{k},s) g^{-1}_{\gamma\alpha} (-k_{\gamma} k_{\delta}) g^{-1}_{\delta\beta} \hat{q}^\dagger(\pmb{k},s) where we have used the identity :math:`\int d\pmb{x} e^{+i(\pmb{k} + \pmb{k}')\cdot \pmb{x}} = \delta(\pmb{k} + \pmb{k}')`. The term inside the :math:`ds` integral is computed for each `s` during the solution of the propagator. Also note that :math:`\pmb{g}^{-1} (-\pmb{k k}) \pmb{g}^{-1}` only needs to be computed once for each call shape :math:`\pmb{h}`. By switching to the Duchs2014 notation described in :class:`model_melt_chi_ab.rst` for a system of many polymers .. math:: \beta \pmb{\Sigma} = \sum_p \frac{C \phi_p}{\alpha_p} \pmb{\Sigma}_p where in this equation the first :math:`\sum_p` (i.e. \\sum_p) is the sum over each polymer `p` and the second :math:`\pmb{\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 :math:`\Sigma` is obtained, :math:`H_V` can be minimized by updating :math:`\pmb{h}` using a forward Euler scheme .. math:: \pmb{h}^{j+1} = \pmb{h}^j - \delta t \lambda_h \left[ \frac{\partial H_V}{\partial \pmb{h}} \right] ^j where :math:`\delta t` is the timestep used to update the :math:`\mu` fields, :math:`\lambda_h` adjusts the rate of the cell updates relative to the fields, and :math:`\frac{\partial H_V}{\partial \pmb{h}} = \beta \pmb{h} \pmb{\Sigma}` (as defined above). Finally, we define the "internal stress operator" .. math:: \tilde{\sigma} \equiv \pmb{h} \Sigma \pmb{h}^T