## Abstract

It is commonly considered that the frustration between the curvature energy and the chain stretching energy plays an important role in the formation of lyotropic liquid crystals in bicontinuous cubic phases. Theoretic and numeric calculations were performed for two extreme cases: parallel surfaces eliminate the variance of the chain length; constant mean curvature surfaces eliminate the variance of the mean curvature. We have implemented a model with Brakke's Surface Evolver which allows a competition between the two variances. The result shows a compromise of the two limiting geometries. With data from real systems, we are able to recover the gyroid–diamond–primitive phase sequence which was observed in experiments.

## 1. Introduction

Liquid crystals are called *lyotropic* if they experience phase transitions by adding or removing a solvent [1]. Typically, the solute in a lyotropic liquid crystal (LLC) is amphiphilic, comprising a hydrophilic head-group and a hydrophobic chain (tail). Upon varying solvent content and/or temperature, LLC can display a rich variety of phases, including the lamellar phase, the hexagonal phase, and the bicontinuous or micellar cubic phases. Of these, the bicontinuous cubic phases are the most complex and interesting. They are observed in cells and organelles, and are believed to play important roles in biological processes [2]. It is generally believed that the bicontinous LLC cubic phases can be described by triply periodic minimal surfaces (TPMSs) [3–6].

A LLC can be treated as a packing of curved amphiphile monolayers and water layers. According to the curvature of the monolayers (measured at the neutral interfaces [7,8]), LLC can be classified into two types [9–11]: with respect to normal vectors pointing from oil to water, the monolayers have positive mean curvature in a *normal* (also known as type I, oil-in-water) system, and negative mean curvature in an *inverse* (also known as type II, water-in-oil) system. In the case of a typical amphiphile–water binary LLC system, an inverse bicontinuous cubic phase consists of a pair of monolayers tail-to-tail (a bilayer) draped over the TPMS [12–14], separating two water channels; while in a normal phase, a water layer following the TPMS separates two channels packed by the hydrophobic tails of the amphiphiles [14].

It is commonly believed that LLC phases arise from the competition between two geometric demands [15–17]: uniform curvature and uniform chain length in the monolayers. Indeed, as the amphiphilic molecules in the monolayers are chemically identical, it is conceivable that they have the same spontaneous curvature [18] and relaxed tail length. However, for most of the LLC phases, the curvature and the chain length cannot be simultaneously uniform [17,19]. Hence the origin of *frustration*, which can be quantified by the variances of mean curvature and of chain length.

Considering both variances at the same time is difficult, hence most investigations treat them separately. The neutral interfaces are either modelled as parallel surfaces of the TPMS (e.g. [11,17,20–23]), eliminating the variance of the chain length; then the Helfrich energy is calculated or measured as the energy of the system. Or, alternatively, the neutral interfaces are assumed to have constant mean curvature (CMC); then the Hooke energy is used as the energy of the system (e.g. [12,13,17,21,24]). The parallel surface model is certainly more straightforward to calculate, but the CMC model makes more sense for the inverse hexagonal and micellar phases [17,19], and is arguably more successful in explaining the phase behaviour involving bicontinuous cubic phases [13].

There has been a constant and strong demand [12,17] for a model that allows the competition between the curvature energy and the chain stretching energy. Here, we present a way of modelling LLC structures in Brakke's Surface Evolver [25] that fulfils this demand. We first demonstrate our method on the inverse hexagonal phase, then apply it to the bicontinuous cubic phases, both inverse and normal, with the geometry of G (gyroid, *Ia*3*d*), D (diamond, *Pn*3*m*) P (primitive, *Im*3*m*) TPMSs. This gives, for the first time, a geometry of LLC that compromises the two extremes: the parallel surfaces and the CMC surfaces.

## 2. Model

Following [13], the surface averaged free energy per hydrophobic chain, denoted by *μ*, comprises two parts:
2.1where
is the Helfrich energy or curvature energy [15,26], and
is the Hooke energy or chain stretching energy. Here, *A* denotes the cross-sectional area of a hydrophobic chain, *L* denotes the chain length, *H* is the interfacial mean curvature, *K* is the interfacial Gaussian curvature, *H*_{0} is the spontaneous mean curvature, *L*_{0} denotes the relaxed chain length, and *κ*_{H}, *κ*_{G} and *κ*_{L} denote the moduli for the energetic contributions from, respectively, the mean curvature, Gaussian curvature, and hydrophobic chain stretching. All these quantities should be measured on or from the *neutral interface*, which is the location within the monolayer where the area is invariant upon isothermal bending [7,8]. The average is over the whole surface, that is .

The contribution from the mean curvature can be rewritten as
2.2where *σ*_{x} = 〈*x*^{2}〉 − 〈*x*〉^{2} denotes the squared variance. Similarly, the contribution from the chain stretch can be decomposed into
2.3The frustration of the system is measured as a weighted sum of the squared variances
which is the quantity to minimize in our model.

This model is certainly a simplification. In particular, we ignore the contribution of the tilt energy, the higher order contributions of the curvatures, as well as the non-local interactions of the monolayers, such as van der Waals interaction and hydration repulsion. Nevertheless, it is hoped that the important physical features are retained.

Brakke's Surface Evolver [25] is a software that minimizes energies of triangulated surfaces subject to constraints and boundary conditions. Apart from the usual surface tension energy, or the area functional, Surface Evolver is able to calculate many other quantities on a surface. Many of these quantities can be included in the energy to minimize. A quantity can be calculated on a particular set of geometric elements (vertices, edges, faces and bodies), allowing the users to control which elements contribute to the total energy. Moreover, each quantity has a modulus to specify its weight in the total energy. The moduli are adjustable in real time on each geometric element.

In particular, Surface Evolver has implemented the Willmore energy (*h* − *h*_{0})^{2} (see [27]) evaluated on vertices, and the Hooke energy (ℓ − ℓ_{0})^{2} evaluated on edges. Here, we use lower case letters to distinguish dimensionless quantities in Surface Evolver from physical quantities. More specifically, our computations in Surface Evolver assume unit lattice parameter (edge length of the conventional cubic cell of the TPMS). If the physical system has lattice parameter *a*, we have the relations *h* = *Ha*, *k* = *Ka*^{2} and ℓ = *L*/*a*. Then, the surface averaged energy per hydrophobic chain takes the following form:
2.4

The parameters *h*_{0} and ℓ_{0} in Surface Evolver are supposed to mean the dimensionless spontaneous mean curvature and relaxed chain length. However, we keep updating these parameters to the average values 〈*h*〉 and 〈ℓ〉, so that the averaged Willmore energy gives the squared variance of mean curvature, and the averaged Hooke energy gives the squared variance of chain length. Hence, the minimized quantity is the frustration
or, for convenience,
where λ = *κ*_{L}*a*^{4}/2*A**κ*_{H}. After sufficiently evolving the surface, it is easy to transform the minimized frustration into the energy using equations (2.2) and (2.3).

Our procedure of modelling bicontinuous LLC cubic phases in Surface Evolver can be outlined in the following five steps:

(1) Prepare a triangulation of the D, P or G surface, or any other surface of interest.

(2) Make two copies of the triangulation, modelling the neutral interfaces. Assign Willmore energy to these copies.

(3) Create edges modelling the hydrocarbon chains, and assign Hooke energy to these edges.

(4) Impose a constraint on the volume

*ϕ*_{n}of the space between the neutral interfaces, in correspondence with the desired water fraction.(5) Evolve the surface towards the minimum of the total energy, while keeping the average values

*h*_{0}= 〈*h*〉 and ℓ_{0}= 〈ℓ〉 updated.

In step (4), the space between the neutral interfaces is the one occupied by the hydrocarbon chains for an inverse system, or the one containing water for a normal system. In step (3), the edges modelling the hydrophobic chains depend on the system and the model. For inverse LLC phases, we are aware of several ways of modelling the chains:

(a) Connect edges between the corresponding vertices on the triangulated neutral interfaces.

(b) Connect edges from the vertices on the neutral interfaces to the corresponding vertices on the TPMS.

(c) Subtending normal vectors from the TPMS until the neutral interfaces (used in [12]).

(d) Subtending normal vectors from the neutral interfaces until the TPMS.

To consider the curvature and the chain length at the same time, we need to keep the chains up to date at every step of surface evolution. The models (c) and (d) require computing the intersection point of a line and a triangulated surface, which is very time consuming. We will use the models (a) and (b). The chains are established once at the beginning, and will be used during the entire computation.

However, only the model (d) forces the chains to be perpendicular to the neutral interfaces. This may raise concerns about tilting (see [28,29]) in other models. In our computations, the tilting is suppressed by the `normal_motion` mode of Surface Evolver, which forces the vertices to move along the normal directions. This practice reduces degrees of freedom and accelerates the computation, and tilting is suppressed as a side effect. However, it does not completely eliminate the tilting. The amount of tilting under the ‘assumption’ of `normal_motion` is measured and presented, but systematic calculation of tilt energy is not within the scope of the current paper.

An edge in the model (a) corresponds to two chains. One advantage of this model is that it is independent of the TPMS: the TPMS is present only for reference, and does not participate in the calculation. Our program is robust in the sense that, if one starts from a periodic CMC surface instead of a TPMS (which the author did, accidentally), the bilayer will correctly evolve to the TPMS. Hence, our program confirms, under the current theory, that bicontinuous LLC cubic phases do follow the geometry of TPMS.

For normal LLC phases, we subtend normal vectors from vertices of the triangulated TPMS to find intersection points on the medial surface, then model the hydrophobic chains by edges from these intersection points to the corresponding vertices on the neutral interfaces. This model follows the idea of [30]. See figure 1 for an illustration of our models.

Readers are free to use other chain models in their own calculations. The technical details are postponed to appendix A.

## 3. Results

### 3.1. Inverse hexagonal phase

In the inverse hexagonal LLC phase, water is contained in the tubes formed by amphiphiles arranged in a hexagonal lattice. Owing to its simplicity, the inverse hexagonal phase has been a standard example (e.g. [17,19]) to illustrate the incompatibility between uniform curvature and uniform chain length, and to demonstrate the calculation of chain frustration within the CMC model.

We also choose the hexagonal phase to explain our method, since the model can be built up in dimension 2. More specifically, the hexagonal lattice graph plays the role of TPMS, and a ‘triangulation’ is nothing but a subdivision of the edges. The neutral interfaces are modelled by cycles within the hexagonal cells. In practice, such a cycle is initially created as a copy of the hexagon. Vertices of the cycles are connected to the corresponding points on the hexagonal lattice graph, modelling the hydrophobic chains of the amphiphiles (model (b)). The area between the cycle and the hexagon corresponds to half of the volume between the neutral interfaces.

Surface Evolver is able to apply squared curvature energy on the cycle, and Hooke energy on the hydrophobic chains. We can choose to minimize either energy alone, or minimize a weighted sum of the two. At each iteration, we update the parameter ℓ_{0} to the average chain length, so that the averaged Hooke energy actually gives the dimensionless squared variance of chain length. Note that the average curvature is a constant here.

In figure 2, we show the result of Surface Evolver under three different situations. Minimizing curvature energy alone results in a circle,^{1} and minimizing chain stretching energy alone results in a star-like shape. If both are present in the total energy, a competition arises. The cycle evolves to a closed curve, which is neither a circle nor a star, but an intermediate shape.

### 3.2. Inverse lyotropic liquid crystal phases

We now repeat the calculation of Shearman *et al.* [12] to verify the feasibility of our method. More specifically, we measure the chain frustration in inverse LLC phases assuming CMC neutral interfaces. The hydrophobic chains are modelled by edges connecting corresponding vertices on the triangulated neutral interfaces; recall that they are copies of the same triangulation. The chain lengths are then half of the edge lengths. The volume fraction *ϕ*_{n} of the space between the neutral interfaces ranges from 0.01 to 0.50. The squared variance *σ*^{2}_{ℓ} of dimensionless chain length is measured as the frustration of the system.

The result is plotted in figure 3*a*. We use 1 − *ϕ*_{n} on the *x*-axis to facilitate the comparison with fig. 5 in [12].^{2} The two plots are very similar: the frustration is the lowest in the G phase, and the highest in the P phase when the water fraction is low, or in the D phase when the water fraction is high. In particular, we also observe a crossover between the P and D phases around *ϕ*_{n} = 0.29 (≈0.25 in [12]). There is only a slight numerical disagreement, which can be explained by the minor differences in our models. In particular, the hydrophobic chains in [12] are perpendicular to the TPMS,^{3} while we allow slight tilting for the chains.

Now that we have confirmed the validity of our algorithm, we turn on the Hooke energy in Surface Evolver, and minimize the weighted sum of the two dimensionless squared variances. We fix the weight of *σ*^{2}_{h} to be 1, and vary the weight of *σ*^{2}_{l}. In other words, we are minimizing *σ*^{2}_{h} + λ*σ*^{2}_{ℓ}, corresponding to the frustration divided by 2*A**κ*_{H}*a*^{−2}. In view of parameter values provided in [13], we choose λ = 10^{3}, 10^{4} and 10^{5} in our computations; see context of figure 5. The resulting surface is neither parallel nor CMC, but a compromise of the two. The measured frustrations, together with the contribution of chain length frustration, is plotted in figure 3.

The mean curvature and the chain length compete to be uniform, and a higher weight is an advantage in this competition, hence the value of *σ*^{2}_{ℓ} is lower comparing with the CMC model. Moreover, the contribution of λ*σ*^{2}_{ℓ} in the total frustration is decreasing as λ increases. This trend is more significant in the D and P phases than in the G phase. As a consequence, the frustration in the G phase eventually surpasses the D and P phases when λ = 10^{5}.

As an example, the evolved P and D neutral interfaces with *ϕ*_{n} = 0.3 and λ = 10^{4} are shown in figure 4. Accompanying histograms show tilt angles less than 5°. The difference of the resulting neutral interfaces from the CMC interfaces with the same volume fraction *ϕ*_{n} is not visible with human eyes (0.5% of the lattice parameter). Hence, we use CloudCompare [32] to measure the deviation from the CMC interfaces, and colour the surface accordingly for visualization. Roughly speaking, comparing to the CMC model, the neutral interfaces tend to move away from the TPMS around the flat points (where Gaussian curvature vanishes) of the TPMS, and towards the TPMS at other places. This is compatible with the observation in the CMC model that the chains are compressed around the flat points of the TPMS, and extended elsewhere; see fig. 8 of [12], or in the electronic supplementary material, figure S1, for a colorful reproduction.

If the water fraction decreases further, the neutral interfaces in the CMC model tend to spheres (positive Gaussian curvature) connected by small necks (negative Gaussian curvature) [33,34]. Eventually, there will be a pinch-off point due to certain physical threshold, e.g. molecular size, and the bicontinuous phases are substituted by micellar phases. For the P phase, the pinch-off is supposed to occur at *ϕ*_{n} = 0.5 [12]. This can be observed in figure 3*a*, as the curve for the P phase becomes very sloped near *ϕ*_{n} = 0.5. However, water fraction in experiment can be much smaller [13]. This inconsistency is a weak point of the CMC model. In our model, the plot curve becomes less sloped, implying that the pinch-off point is moving towards lower water fraction. Indeed, we see in figure 4 that the competition causes expansions in the necks, therefore delays the pinching off. Hence, our model is able to cover highly dehydrated LLCs observed in experiment, which cannot be explained by the CMC model [13,20].

Calculations shown in figure 3 assume constant modulus, thus do not correspond to any experimental systems. To connect to real systems, we need to know the lattice parameter *a*. Recall that the surface averaged energy per hydrophobic chain divided by *κ*_{H}, expressed with dimensionless quantities, is
The total Gaussian curvature is a topological constant 2*π**χ*, where *χ* is the Euler characteristic of the TPMS in the conventional cubic unit. The value of *χ* is −4 for the P surface, −2 for the D surface and −8 for the G surface. An expression of *a* in terms of *ϕ*_{n} is obtained in [21]; typical ranges for other parameters are available in [13]. We take the values *A* = 33 Å^{2}, *κ*_{G}/*κ*_{H} = −0.75, *κ*_{L}/*κ*_{H} = 0.00035 Å^{−2}, *H*_{0} = 1/62.8 Å^{−1}. As for the relaxed chain length, we take *for the moment* *L*_{0} = 8.8 Å, which is the measured distance from the chain ends to the neutral interface in the G phase [35] of 1-monoolein/water system. These allow us to plot *μ*/*κ*_{H} in figure 5 as a function of 1 − *ϕ*_{n}. By taking 5° as the typical tilt angle, and *κ*_{θ}/*κ*_{H} = 10^{18} for the tilting modulus [28,29], we can estimate that the tilt energy is smaller by one order of magnitude.

The phases of lowest energy give a G–D–P sequence with increasing water fraction, which is a widely observed phase sequence in experiments [36,37]. Meanwhile, disagreements with experimental data [38] are also visible. In particular, the P phase is not present in the monoolein/water system; and the phase transitions seem to occur at a higher water fraction compared with experiments.

We point out the major difficulty at the current stage. Our program can achieve a very high precision (*σ*^{2}_{H} <10^{−26} for CMC surfaces). However, the relation between *a* and *ϕ*_{n}, as well as some parameter values we use in the calculation above, are actually derived under the CMC model. Despite the closeness in some statistics, our model differs significantly from the CMC model at a low water fraction. At *ϕ*_{n} = 0.5, the deviations of the evolved neutral interfaces from CMC interfaces become comparable with the chain length (electronic supplementary material, figure S2). In a complicated system as the one we are looking at, minor energy differences could lead to different phase behavior. Better estimates of the parameters would be crucial for a better prediction.

### 3.3. Normal lyotropic liquid crystal phases

LLC in normal bicontinuous cubic phases seem to be less common. The authors are not aware of any report other than normal bicontinuous G phase [39,40]. Meanwhile, in mesophased silica system, which can be seen as LLC systems with water replaced by silica, various normal phases have been reported, including the bicontinuous G [41–43], P [44] and D [45] phases.

Despite our limited understanding, we repeat the same calculations to the normal bicontinuous LLC cubic phases, with the volume fraction *ϕ*_{n} of the space between the neutral interfaces (containing water) ranging from 0.01 to 0.50. Following Schröder-Turk *et al.* [14,30], the hydrophobic chains are modelled by edges from vertices on the triangulated neutral interfaces along the normal vectors to the nearest point on the medial surface of the TPMS.

Because of the way that Surface Evolver calculates the volume, the measurements near *ϕ*_{n} = 0 are not physically meaningful; see appendix A. The tilt angles are much larger than in the inverse phases, hence the tilt energy could be quite significant. Hence, we do not attempt to connect our calculation with real systems. For these reasons, and also due to lack of CPU power, we only compare the CMC model with the case λ = 10^{4}.

The measured frustrations for the P, D and G phases, as well as the contribution of chain stretching energy when λ = 10^{4}, are plotted in figure 6. In both CMC model and our model, the frustration is the highest in the P phase and the lowest in the G phase. This could explain the difficulty of observing the normal bicontinuous P and D phases in experiment. As in the inverse phases, we see that the value of *σ*^{2}_{ℓ}, as well as the contribution of λ*σ*^{2}_{ℓ}, declines as λ increases. The effect is again more significant in the D and P phases than in the G phase.

Evolved D and P bilayers with *ϕ*_{n} = 0.3 are shown in figure 7, to compare with figure 4. Again, the neutral interfaces tend to move away from the TPMS around the flat points, and towards the TPMS at other places. Indeed, the flat points of the TPMS are most distant from the medial surface [30]. As a consequence, the pinching-off necks get expanded. Comparing with the inverse phases, the deviation is much larger (greater than 3% of the lattice parameter).

## 4. Conclusion

We have demonstrated that modelling the competition between the curvature energy and the chain stretching energy in LLC systems is possible with Surface Evolver. More specifically, we applied Willmore energy to triangulated surfaces modelling the neutral interfaces, and Hooke energy to edges modelling the hydrophobic chains. A weighted sum of the two is used as the total frustration of the system, and is minimized by Surface Evolver. The resulting surfaces are neither parallel surfaces nor CMC surfaces, but a compromise of the two. A detailed procedure is described and tested for LLCs in bicontinuous cubic phases.

We compare our results on inverse cubic LLC phases to similar calculations with the CMC model [12]. The squared variance of chain length, as well as its contribution of the total frustration, is reduced as the modulus of chain stretching energy increases. This effect is more significant in the P and D phases than in the G phase. The frustration of the inverse G phase eventually exceeds the inverse D and P phases. A closer look reveals that the pinching-off necks in the CMC model, which was the main obstacle to achieve a high dehydration [13,20], get expanded in our model. These observations prove the validity and usefulness of our model. When real data are applied, our calculation yields a G–D–P phase sequence with increasing water fraction, which has been observed in experiments. We also performed same calculations on normal cubic LLC phases. Similar observations are made, but the deviation from the CMC model is much more significant than in the inverse phases.

## Competing interests

We declare we have no competing interests.

## Funding

No funding has been received for this article.

## Acknowledgments

The authors appreciate discussions with Karsten Große-Brauckmann and Ken Brakke, and thank Gerd Schröder-Turk and Lu Han for feedbacks. C.J. acknowledges the support of Corinna C. Maaß and Stephan Herminghaus. We are also grateful to the anonymous reviewers for their valuable comments on a preliminary version of the manuscript.

## Appendix A. Technical details

In the opinion of the authors, the modelling capacity of Surface Evolver is underestimated. One purpose of this paper is to present a non-trivial application of Surface Evolver in the study of interfaces. Our program is flexible and easy to adapt for different scenarios. In this appendix, we provide details of our calculation on bicontinuous LLC cubic phases, as a reference for readers who are interested in carrying out similar computations. The datafiles, source codes, dump of evolved surfaces, as well as the measured data, are publicly available at https://github.com/Dr-How/CubicLLC.

**A.1. Surface preparation**

The cubical periodic minimal surfaces considered in this paper, namely Schwarz G, D and P surfaces, can be generated from a small ‘fundamental patch’ (Flächenstück) by Euclidean motions like translations, reflections, rotations and rotoreflections. However, only translations and reflections apply to the accompanied neutral interfaces that are of physical interest. Reflectional symmetries are related to free boundary conditions, and translational symmetries to periodic boundary conditions. Rotational symmetries are often related to fixed boundary conditions, but the neutral interfaces are not subject to any fixed boundary condition. Hence, we only consider the fundamental unit of the translational and reflectional symmetries.

For the P and D surfaces, the `.fe` datafiles are available on Brakke's website (http://facstaff.susqu.edu/brakke/evolver/examples/periodic/periodic.html). The P patch in the datafile is directly usable. We have to apply a rotation to the D patch to obtain the desired unit. For the G surface, we use a 96-facets datafile kindly provided by Große-Brauckmann [46]. The G surface has no reflectional symmetry, hence the fundamental unit is a translational unit cell. In Surface Evolver, the `torus` model is used for the G surface to impose periodic boundary condition.

After loading the datafile into Surface Evolver, we use the command `r` (refine) to subdivide each triangle into four and, after each refinement, evolve the surface towards the minimum of the area functional. A translational unit cell is not a minimizer of the area functional, unless a volume constraint is imposed [47], which we did for the G surface. Apart from the command `g` that performs one iteration of gradient descent, the command `hessian_seek` applies Newton–Raphson method to the gradient of energy. It is more efficient, yet safer than the more radical `hessian` command.

As in [13], the P patch is refined four times, yielding 1024 faces per patch; and the D patch is refined five times, yielding 2048 faces per patch. Owing to lack of CPU power, we only refine the G surface three times, yielding 6144 faces per unit cell.^{4} Surface Evolver provides commands to eliminate elongated triangles (`K`) and extreme edges (`l` and `t`), as well as commands for vertex averaging (`V`) and equitriangulation (`u`). They are frequently employed to keep the triangulation in good shape.

Surface Evolver is very good at preparing triangulations. But if an exact formula is known for the surface of interest, one could also use a mesh generator to produce a triangulation of high quality; see for instance [30].

**A.2. Setting up the bilayer**

Now that we are in possession of a triangulation of high quality, we make two copies of it to model the neutral interfaces. The commands `new_vertex`, `new_edge` and `new_facet` are useful in this step. We also use element attributes to remember the correspondence between the elements in the original triangulations and in the copies.

To model an inverse phase, we only need to create edges between the corresponding vertices in the copied triangulations (model (a) as illustrated on figure 1*a*). The edge length is thus twice the length of the chain we intend to model. This practice aims to eliminate the interference of the original TPMS in our calculation. After creation of the chains, the TPMS is fixed, and serves only as a reference. This procedure also works on unbalanced TPMS, such as the I-WP surface.

To model a normal LLC phase, we need to find, for each vertex in the triangulated TPMS, the nearest point on the medial surface along the normal vector. Efficient algorithms for this purpose based on the Voronoi diagram was proposed by Amenta *et al.* [48] and Schröder *et al.* [30]. In particular, [30] contains an extremely detailed description of the medial surfaces of the D, P and G surfaces, accompanied by fantastic figures. The medial surface of the P surface is very simple: they are contained in the reflection planes. Hence, we simply attach one end of the edge to a vertex *v* on the P surface, extend the edge along the normal vector at *v* (available in the `vertex_normal` attribute), and attach the other end to the first point we encounter on the reflection plane. The same practice also yields a good approximation for the medial surface of the D surface.

For the G surface, we implement with Scipy a slightly simplified version of the algorithm in [48]. For a vertex of the triangulated G surface, we choose from its Voronoi cell the vertex with the largest inner product with its normal vector. This is a good approximation of the ‘pole’ in [48], which is defined as the furthest vertex of the Voronoi cell. Scipy computes Voronoi diagrams using the Qhull library [49].

Setting up for the inverse phases is automated with the script language of Surface Evolver. Automation for the normal phases is possible with exterior programs.

**A.3. Minimization and measurement**

In Surface Evolver, the Hooke energy (ℓ − ℓ_{0})^{2} is implemented as `hooke_energy`. The quantity `star_perp_sq_mean_curvature` is a robust implementation of the Willmore energy (*h* − *h*_{0})^{2}; see [27]. The Willmore energy is internally weighted by the effective area of each vertex, which is 1/3 of the total area of facets containing the vertex. The parameters *h*_{0} and ℓ_{0} are meant to be the spontaneous mean curvature and the relaxed length. We keep updating them to the average values 〈*h*〉 and 〈ℓ〉 so that the energies measure the squared variances. The surface tension energy is turned off (`set facet tension 0`).

The moduli of the energies are also used for taking averages and normalizations. Let *κ*_{L} and *κ*_{H} be the experimental moduli for the chain stretching and mean curvature energy, respectively. Then in Surface Evolver, the modulus for the chain stretching energy is *k*_{L}*a*^{2}*A*_{v}/*a*^{2}_{0}*A* at vertex *v*, and the modulus for the curvature energy is *k*_{H}*a*^{2}_{0}/*Aa*^{2} for the dimensionless squared variance. Here, *A* is the total area; *A*_{v} is the effective local area of the vertex *v*; *a* is the experimental lattice parameter; *a*_{0} is the lattice parameter in Surface Evolver model, which is 2 for the P and D surfaces, and 8 for the G surface. Recall that the Willmore energy is internally weighted by effective area, hence a factor 1/*A* suffices for normalization.

Let *ϕ*_{n} be the desired fraction of the space between the neutral interfaces. We impose a volume constraint of value *ϕ*_{n}*V* to the body bounded by the neutral interfaces. Here, *V* is the volume of the fundamental unit, which is for P, for D and 256 for G surface. For normal phases, the neutral interfaces may intersect when *ϕ*_{n} is close to 0. This is physically impossible, but not a numerical mistake. Surface Evolver computes volume by vector integration on the bounding faces. In our situation, the integrations on the two neutral interfaces have opposite signs, and they cancel to the imposed small volume. As *ϕ*_{n} increases, the neutral interfaces will be separated very soon.

In practice, the volume fraction *ϕ*_{n} is changed gradually in steps of 0.01. The small changes not only ensure a precise measurement, but also avoid brutal changes in the surface that could destroy the surface. After each step, the surface is evolved (by commands `g` and `hessian_seek`) towards the minimum of the total energy. To measure the dimensionless squared variances, the parameters *h*_{0} and ℓ_{0} are updated after each iteration to the area weighted average values 〈*h*〉 and 〈ℓ〉. After the surface stabilizes, the values of the energies are printed on screen by the command `Q`, or output to an exterior file for future record.

To measure the chain stretching frustration in the CMC setting, we could simply change the quantity `hooke_energy` from `energy` to `info_only`. Then, the Hooke energy is calculated, but excluded from the energy to minimize. For the G phase, however, lack of free boundary condition allows the surface to shift away in the absence of `hooke_energy`. This would increase the variance of chain length. To avoid this, it is recommended to start with non-zero modulus of `hooke_energy`, change the modulus to 0 after several iterations, and change it back for measurements after obtaining a CMC. In practice, we are usually able to reduce the squared variance of mean curvature down to the order of 10^{−26} or lower, hence the obtained surface is indeed of CMC.

**A.4. Performance**

The program is very robust. The author once accidentally applied the program to a CMC gyroid, then the neutral interfaces evolved away from the original position, and correctly stabilized beside the minimal gyroid. In practice, this means that one only needs a general idea of the surface to carry out our method. Meanwhile, brutal changes are still to be avoided: on higher dimensional problems, gradient descent method often suffers from saddle points and local minimums. Hence, we recommend incremental change of *ϕ*_{n} at small steps.

Tilt of the chains is allowed, but is suppressed by the `normal_motion` and `hessian_normal` modes, which forces the vertices to move along the normal directions of the surface. Tilting is very minor for inverse phases, but more significant in normal phases. We do not pay much attention to the tilt energy, but the amount of tilting is measured and presented for reference. With some more work, it is possible to include the tilt into the energy to minimize. This is, however, not our current focus.

Distortion of the triangulation is the main source of inaccuracy. Elongated triangles and edges of extreme lengths would cause problems to Surface Evolver. An initial triangulation of high quality is recommended for this reason. But the problem is unavoidable when the bilayer becomes very thick. In our calculations, the measurement suffers from slowness and imprecision as *ϕ*_{n} approaches 0.50. Surface Evolver provides commands to modify the triangulation and improve the quality. But most of these commands involve vertex deletion, which would destroy elastic edges. Our options are limited to two commands: `V` for averaging vertices, and `u` for equitriangulation.

The measurement can be automated, but human supervision is recommended to achieve a good precision. The `dump` command saves intermediate status of the surface, allowing the user to check for problems and perform additional iterations to improve the precision. A full measurement of G phase would take a few hours on a personal laptop.

## Footnotes

One contribution of 17 to a theme issue ‘Growth and function of complex forms in biological tissue and synthetic self-assembly’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3776552.

↵1 This is related to a quite challenging mathematical problem [31].

↵2 The ‘

*ϕ*_{w}’ in [12] should be 1 −*ϕ*_{n}. The relation between*ϕ*_{n}and*ϕ*_{w}is*ϕ*_{n}= (1 −*ϕ*_{w})*ν*_{n}/*ν*, where*ν*is the molecular volume and*ν*_{n}is the volume between the neutral interface to the chain ends.↵3 To be rigorous, the chains should be perpendicular to the CMCs.

↵4 Shearman

*et al.*also claim to have refined three times in [12], but report 24 576 facets. Their program provided as the electronic supplementary material shows four refinements.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.