Approximate Modeling of Spherical Membrane

Spherical symmetry is ubiquitous in nature. It's therefore unfortunate that spherical system simulations are so hard, and require complete spheres with millions of interacting particles. Here we introduce an approach to model spherical systems, using revised periodic boundary conditions adapted to spherical symmetry. Method reduces computational costs by orders of magnitude, and is applicable for both solid and liquid membranes, provided the curvature is sufficiently small. We demonstrate the method by calculating the bending and Gaussian curvature moduli of single- and multi-layer graphene. Method works with any interaction (ab initio, classical interactions), with any approach (molecular dynamics, Monte Carlo), and with applications ranging from science to engineering, from liquid to solid membranes, from bubbles to balloons.


I. INTRODUCTION TO MODELING APPROACH
The problem in simulating spherical symmetry is topological: you cannot build a perfect sphere from identical blocks. The absence of such a building block has enforced expensive simulations with complete spheresthough usually spherical simulations are simply avoided. Overwhelming dilemmas like this are often considered so fundamental and frustrating that they restrain all attempts to seek for a practical solution.
Anyhow, avoiding spherical systems in our world is hard. Spherical shells surround us in a variety of forms: in balloons, in cell membranes inside our bodies, in bubbles in the sea, or in Earth's crust. The interaction of nanoparticles with cell membranes, for instance, is a topical question. 1 Since cell membranes' curvature moduli determine the very forms of red blood cells, for example, one can see why simulations should incorporate curvature effects. 2,3 Another timely example is the foam of spherical bubbles in the sea, the bursting of which may play an important role on the so-called sea spray that produces spherical aerosols into the atmosphere. 4,5 Although liquid membranes are more abundant in nature, also man-made solid membranes have spherical symmetries, at least locally. Examples are fullerenes, 6 nanoballoons, 7 and especially graphene that contains intrinsic ripples even when suspended freely. [8][9][10] Curvature moduli of graphene are intimately related to these ripples, whether they are intrinsic or not, 11,12 and in a broad sense to elastic behavior of all honeycomb carbon, among graphene nanoribbons, 10 multilayer graphene, 13 and carbon nanotubes. 14,15 Conventionally spherical systems are treated in three ways. The first way is to simulate the system as a whole. Needless to say, this is expensive and often impossible. 3 The second way is, should the system have some welldefined point-group symmetries, to use those symmetries for reducing computational costs. Most established codes have the ability to benefit from such symmetries; this has long been a standard procedure with molecules and clusters. 16 Because the symmetry is exact, however, neither the curvature nor other geometrical parameters can * Corresponding author be changed flexibly. The third way is to ignore curvature altogether and to use periodic boundary conditions (PBC) to simulate an infinitely large, flat membrane. Unfortunately, in nanoscience many systems fall between these two extremes: systems with huge number of particles, having no overall symmetry, but prominent curvature effects. At the moment a practical way to simulate such systems does not exist.
The periodic boundary conditions have been adapted, however, also to symmetries beyond translation. The first ideas came along chiral carbon nanotubes, 17,18 and those ideas have been used ever since; for reviews look at Refs. 19 and 20. An important extension to general symmetries with exact treatment was done in Ref. 21, which has enabled more flexibility. [22][23][24] . Later, Ref. 25 introduced revised periodic boundary conditions (RPBC), a simple formalism for general material distortions; this is the approach we shall use here, and it's illustrating to review it briefly.
In RPBC, the usual translation operations are replaced by general symmetry operations S n that, in a quantummechanical language, leave the electronic potential invariant, orD The operation S n is a succession of an abelian group of operations S i , that is S n = S n1 1 S n2 2 · · · . Then, by imposing periodicity (S Mi i = 1, M i integer), one finds that the Hamiltonian eigenstates ψ aκ (r) at r and at r = S −n r differ only by a phase factor, D(S n )ψ aκ (r) = ψ aκ (S −n r) = exp(iκ · n)ψ aκ (r), (2) with inverse operation S −n , band index a, and the reciprocal lattice vector κ. Eq.(2) infers the familiar result: a single simulation cell-whatever its shape-is enough to describe the extended system as a whole. Revised PBC is hence similar to conventional PBC and differs only in the definitions of the symmetry operations. There are no other fundamental differences. As an illustrative example, the total energy with a classical pair potential is where N is the particle count and n runs over operations where particle i at R i still interacts with the periodic arXiv:1010.0067v2 [cond-mat.mtrl-sci] 23 Nov 2010 In this paper we use the RPBC, reviewed above, in an approximate way to introduce a trick for modeling spherical membranes. Adapting RPBC for spherical systems enables simulations with orders-of-magnitude reductions in computational costs. We shall apply the method to calculate graphene's mean and Gaussian curvature moduli, but first we proceed to discuss symmetry operations and their character.

II. SPHERICITY AS AN APPROXIMATE SYMMETRY
Consider the square cone in Fig. 1, regard the grid as fixed in space, and concentrate on the shaded region. If we rotate all particles an angle δθ 1 around y-axis, or an angle δθ 2 around x-axis, the geometry within the shaded region will remain approximately intact. This means that rotations S n1 1 r = R(n 1 δθ 1 c 1 )r and S n2 2 r = R(n 2 δθ 2 c 2 )r (with c 1 =, c 2 =î and operation R(c) as |c|-radian rotation aroundĉ) leave the electronic potential V (r) invariant near the shaded region: S 1 and S 2 are symmetry operations as far as the shaded region and its vicinity is concerned. Two rotations around different axes do not commute in general, but if the rotation angles δθ i are small Hence also the combined operation is an approximate symmetry operation, provided that n i are small enough. Eq. (4) is basically all we need to fully employ the RPBC of Ref. 25; we are, in principle, ready to go and to simulate any spherical membrane.

III. FEATURES DUE TO APPROXIMATION
In practice, however, the approximate character of S n raises questions that deserve some elaboration. First, as already mentioned, the formalism assumes periodic boundary conditions (S Mi i = 1) which may seem questionable. Here we remind that similar PBCs are used also in regular bulk, with all three dimensions periodic in an intertwined fashion. (In two dimensions PBC represents topologically a toroid.) The bottom line is that periodicity is not a physical reality but a mere mathematical trick that works, and enables the application of revised Bloch's theorem in the first place. 27,28 The integers M i are connected rather to κ-point sampling than to physical reality.
Second, revised PBC does not need the "unit cell" concept. However, we shall call the square cone in Fig. 1, extending from the origin to infinity and enclosing the shaded region, a unit or simulation cell because the concept is familiar and convenient in discussion. Otherwise, the mere expression for S n in Eq.(4) is enough to determine everything in the simulation.
Third, the claim is not to simulate a complete sphere, but rather to view the curvature as a local property. The particles in the simulation cell see the closest environment curved-and only this is important. The simulation cell is the only cell we model, and distances and angles measured only from the simulation cell are meaningful. For example, the vicinity of particle at r in Fig.1 exhibits curvature in bond angles and distances if one looks at particle's own periodic images at S 2 r and S −1 2 r. Symmetry operations S n that have n i large enough to rotate large angles (n i ∼ (π/2)/δθ i ) should be excluded because the non-commutativity of S i 's would otherwise become significant.
Fourth, the radius or curvature R in Fig. 1 is not a parameter in the simulation; radially particles can migrate wherever interactions drive them. The spherical form is only forced by the choice of symmetry operations and the parameters δθ 1 and δθ 2 , and since the symmetry is discrete, the system needs to be neither continuously nor smoothly spherical.
Fifth, a natural limitation is to have enough empty space near the origin to avoid too close encounters between the particles. 29 Membrane can be thick.

IV. APPLYING THE METHOD
The validity of the method depends on the system and its interactions. As a principal rule, the radius of curvature R should be much larger than the interaction ranges between the particles. If ranges are larger than the system size, especially if those interactions control morphology, one does better to model the complete system. The Coulomb interactions can play a role locally, within small length scales (size of the unit cell at most), but the longranged Coulomb interaction requires special care, perhaps some refinements (the unit cell better be neutral). 30 Quantitative error due to the non-commutativity of S i 's can be estimated by first using the right-hand side of Eq.(4) as S n , and then using the left-hand side of Eq.(4) as S n (changing the ordering of S n1 1 S n2 2 ), and comparing the resulting energies.
Because liquid lacks long-range order, the method suits particularly well for liquid membranes, such as lipid bilayers. Their energetics can be described by the Helfrich Hamiltonian that gives membrane's elastic energy per unit area as 31 Here κ is the mean curvature modulus (don't confuse with a κ-point), κ is the Gaussian curvature modulus, and R 1 and R 2 are the principal radii of curvature. The liquid membrane doesn't need to be free-standing, however, because also solid support can be incorporated, either by external force fields or by fixed atoms. External radial forces can be also used for pressurization, mimicking the embedding of membrane in gaseous or liquid environments. For solid membranes the situation is more complicated, because energy will come also from the internal strain E s . If a flat, round sheet of radius ρ is wrapped into a spherical segment, the energy will be E s ∼ Ehπρ 6 /108R 4 , 32 where E is the Young's modulus of the material, h is the membrane thickness, and R is the radius of curvature; meanwhile the curvature-related energy is E c = g · πρ 2 . Hence, for a reasonable modeling of solid membranes using Eq.(5), we need to have E s E c , or which suggests a minimum radius of curvature R min for a given unit cell area. If this geometrical and materialdependent criterion should be violated, the simulation would be dominated by non-local stress fields. Since the method does not properly describe these fields, the treatment would become ill-defined. The above problem is present when sphericity is forced on originally flat sheet. But defects, for example, can induce spontaneous curvature in solid membranes in which case R min can be smaller. The method provides a new tool to investigate phenomena such as rippling due to adsorption-induced pinching of the membrane. 11 This method does not directly compete with any existing method, but instead it provides possibilities to do something new.

V. EXAMPLE: SPHERICAL GRAPHENE
The spherical symmetry was implemented in the density-functional tight-binding code hotbit. 33,34 The RPBC implementation has a negligible computational overhead as compared to translational symmetry, 25 and can be implemented just by a few lines of new code in any existing RPBC implementation. The code source is open and stands for inspection.
In this section we use the hotbit implementation to present one practical example. We calculate the curvature moduli of graphene, motivated by their relevance to present-day engineering with carbon nanostructures. For a sphere the radii of curvature are R 1 = R 2 = R, and Eq.(5) gives g = (2κ + κ)/R 2 ; for a cylinder R 1 = R, R 2 = ∞, and g = κ/(2R 2 ). Hence, by calculating the elastic energies for a cylinder and a sphere and varying δθ i 's (hence varying R) we obtain both κ and κ directly.
Prior to simulating spherical graphene, we first calculated the mean curvature modulus of graphene, also applying revised PBC. Only now the symmetry operations, in a cylinder-like setup, were a rotation around z-axis (S 1 ) and translation in z-direction (S 2 ) with a 4atom unit cell (like a nanotube with enormous diameter); we won't discuss the cylinder setup further here. 18 The resulting cohesive energy depends on R quantitatively like R −2 , as Eq.(5) suggests, and the fitted value for κ = 1.61 eV (4.22 eVÅ 2 /atom) agrees with a densityfunctional reference value (1.5 eV) 35 albeit is larger than an experimental reference value (1.2 eV). 36 Returning to spherical graphene, Figs. 2a and 2b show the two-atom unit cell of graphene. Unlike in Fig.1, the unit cell is skewed with c 1 = and c 2 = cos(5π/6)î + sin(5π/6). The geometry was optimized with given δθ i 's, which were taken as δθ i = 2.5Å/R when we wanted to investigate a radius of curvature that roughly equals R . 37 All the radii of curvature we report, anyhow, are the optimized R (R ≈ R because curvature changes bond distances only slightly). In practice we found that structure optimizations require convergence criteria tighter than with translational cells, due to geometrical effects from small δθ i . 38 In quantum simulations κ-points can be freely sampled (κ i ∈ [−π, π]) because PBC is an ap- proximation, just as with conventional Bloch's theorem; we used a 50 × 50 κ-point mesh. Fig.2c shows our main result, graphene's cohesive energy as a function curvature-and represents the showcase of the new physics this method can unearth. Energy behaves clearly like R −2 , as suggested by Eq.(5). The energy penalty 6.6 eVÅ 2 R −2 /atom, combined with previously calculated κ, yields κ = −0.70 eV; we could not find this number in the literature. This result confirms graphene's beautiful elastic behavior up to high curvature-also for spherical distortion. 35 We did consistency checks for the graphene sphere calculations, three listed next. As a first check, when we investigate Eq.(6) with graphene parameters, we get ρ √ 6Å · R. For a graphene unit cell ρ ∼ 1Å (lattice constant 2.5Å), and the consequent criterion R 0.2Å is easily fulfilled. We obtained the same κ with N = 8 and N = 32 atom unit cells, even though larger N increases R min (Eq.(6) and ρ 2 ∝ N infer R min ∝ N ). Thus, the area is small enough to be stress-free, and the simulation is indeed dominated by curvature energy alone. We were able to perform controlled calculations down to radii R min ∼ 10Å or δθ max ∼ 15 • . As a second check, we estimated quantitative error in energy due to the noncommutativity of the two rotations (inset in Fig.2c), as suggested above, but found the error fairly small. As a third check, we implemented symmetry also with a negative Gaussian curvature R 1 = −R 2 = R, for which g = −κ/R 2 directly, and got an independent confirmation for κ; we won't attempt to describe structures with negative Gaussian curvature here. Finally, since there is no charge transfer, the long-range Coulomb interactions are no issue.
Closer inspection of geometry revealed that curvature increased bond distances as d nn = 1.417Å+0.135Å 3 /R 2 , due to the weakening of in-plane σ-bonds, and hereby decreasing the effective nearest-neighbor tight-binding hopping parameter as t eff = t gr − 4.8 eVÅ 2 /R 2 (t gr ≈ 2.7 eV). For a detailed discussion of curvature-induced effects on graphene, we recommend Refs. 39, 40 and 41.
For completeness we calculated κ and κ for bi-and trilayer graphene as well, and summarize the results in Table I. Assuming a constant layer separation of h = 3.4Å analytical expressions for the curvature moduli of multi-layer graphene come as κ n =nκ 1 + Eh 3 (n 3 − n)/12 κ n =nκ 1 − Eh 3 (n 3 − n)/12, where n is the number of layers and E is Young's modulus. The simulated and analytical numbers have a fair agreement. Table reveals how strikingly smaller the moduli are for graphene monolayer, a true oddity among solid elastic sheets, as noted already in Ref. 42.

VI. CONCLUDING REMARKS
We have introduced a simple and practical method to simulate spherical systems using revised PBC. Although the method is approximate, it is applicable precisely to systems so hard to handle: large systems with prominent curvature effects. Since the method works with schemes from ab initio electronic structures and classical potentials to coarse-grained and finite element modeling, and has a wide range of applicability, we encourage any additional implementations.
Admittedly, it may take some time to digest the approximate nature of the method. The role of symmetries in materials modeling is usually taken as clear-cut, solid, and untouchable: it either is or is not. In this paper we have, however, created and entered a new gray area in symmetry usage; we are unaware of symmetry being treated in this type of approximate fashion before. For this reason, when using approximate spherical symmetry-or other approximate symmetries in futurewe urge to examine modeled systems carefully and get assured of method's validity; the best guide on this way is common sense.