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τ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 ee1, ee2, ee3, we define the cell shape tensor

hh=(ee1ee2ee3)

The metric tensor is defined as

gghhThh

To find the optimum box size, we seek to minimize HVhh where HVH/V is the intensive energy density. This derivative is related to the internal stress operator Vσσint=HhhhhT. For a general multi-chain Hamiltonian, HHlocpCV~ϕpαplnQp. The stress operator is thus

Vσσint=HlochhhhTpCVϕpαpQpQphhhhT

For most of the polymer models described in GHF (i.e. without Coulomb interactions), Hlochh=0 and the calculation of Vσσint involves the calculation of Qphh. To derive this expression, one first shows that Qphh=2hhQpgg. Recalling that Qp=dxxqp(xx,αp)=1Vdrrqp(rr,αp) where qp satisfies the general PDE sq(xx,s)=Lq(xx,s), one can derive

Qpgg=dxx0αpqp(xx,αps)Lggqp(xx,s).

For a continuous Gaussian chain model, L=iw(xx)=gg1:xxiw(xx) and so

Qpgg=dxx0αpqp(xx,αps)hh1ThhTqp(xx,s).

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

Vσσint=HlochhhhT+p2CVϕpαpQp1Vdrr0αpdsqp(rr,αps)Tqp(rr,s)

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

Vσσint=p2CVϕpαpQpdrr0αpdskkqp(kk,αps)(kkkkT)qp(kk,s)

where we have now set Hlochh=0.

For smeared models where Qp[Γw], an addition term emerges since Γ depends on gg. For smeared models the stress is

Vσσint=p2CVϕpαpQpdrr0αpdskkqp(kk,αps)(kkkkT)qp(kk,s)+pCVi=1Saikkρ^i(kk)Γ^(kk)w^i(kk)(kkkkT)

where the sum over i runs over each species, ai is the smearing range, and ρi is the the density operator,

ρi(xx)=ϕpαpQpsidsq(xx,αps)q(xx,s).

Once σint is obtained, HV can be minimized by updating hh using a forward Euler scheme

hhj+1=hhjδtλh[HVhh]j

where δt is the timestep used to update the μ fields, λh adjusts the rate of the cell updates relative to the fields, and HVhh=σinthhT (as defined above).

LEGACY DOCUMENTATION

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

The only term in HV that depends on hh is the logQ term. For these models

HVhh=βhhΣΣ

where Σ is a symmetric tensor that for molecule p whose elements Σp,αβ are given by

βΣp,αβ=np3VQdxx0Npds[b(s)]2gγα1q(xx,s)2xγxδq(xx,s)gδβ1

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

Note

The prefactor in Barrat2005 Eq 36 is nV2Rg,02Q 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 Σ is incorrect in the code. However, since we’re usually interested in where Σ=0, it shouldn’t matter (at least for the moment).

Since Σp involves q(xx,s) and q(xx,s) it is straightforward to calculate Σp while solving the modified diffusion equation for the propagators.

Following a computation trick described by Kris Delaney, Σ can efficiently be substituting the definition of the fourier transform of q(rr,s)=1Vkkq^(kk,s)eikkxx. Also recalling that 2xγxδkγkδ in fourier space,

βΣp,αβ=np3VQp0Npds[b(s)]21V2kkq^(kk,s)gγα1(kγkδ)gδβ1q^(kk,s)

where we have used the identity dxxe+i(kk+kk)xx=δ(kk+kk). The term inside the ds integral is computed for each s during the solution of the propagator. Also note that gg1(kkkk)gg1 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

βΣΣ=pCϕpαpΣΣp

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

Once Σ is obtained, HV can be minimized by updating hh using a forward Euler scheme

hhj+1=hhjδtλh[HVhh]j

where δt is the timestep used to update the μ fields, λh adjusts the rate of the cell updates relative to the fields, and HVhh=βhhΣΣ (as defined above).

Finally, we define the “internal stress operator”

σ~hhΣhhT