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