Efficient approach for simulating distorted materials

The operation principles of nanoscale devices are based upon both electronic and mechanical properties of materials. Because these properties can be coupled, they need to be investigated simultaneously. At this moment, however, the electronic structure calculations with custom-made long-range mechanical distortions are impossible, or expensive at best. Here we present a unified formalism to solve exactly the electronic structures of nanomaterials with versatile distortions. We illustrate the formalism by investigating twisted armchair graphene nanoribbons with the least possible number of atoms. Apart from enabling versatile material distortions, the formalism is capable of reducing computational costs orders of magnitude in various areas of science and engineering.

The formalism is obtained by revising Bloch's theorem, and the derivation is straightforward. Consider electrons in a potential V (r) that remains invariant in symmetry operations S n , D(S n )V (r) = V (S −n r) = V (r). (1) The operation S n , with inverse S −n , is a succession of n i times operation S i for all i = 1, 2, . . ., that is S n = S n1 1 S n2 2 · · · with n = (n 1 , n 2 , · · · ). S i can be any symmetry operation, such as translation, rotation, reflection, inversion, joined translation+rotation, or joined translation+reflection, to mention six, and they should form an abelian group.
Let us give a couple of familiar examples. With bulk S i 's are three translations; for a benzene ring (C 6 H 6 ) S 1 could be a two-, three-, or six-fold rotation around the symmetry axis; for an achiral carbon nanotube S 1 could be a translation along the symmetry axis and S 2 could be, say, an M -fold rotation around the symmetry axis; for polyethene ([-CH 2 CH 2 -] n ) S 1 could be a translation across one CH 2 unit followed by a reflection (normally S 1 would be a translation across the whole CH 2 CH 2 unit). Now, returning to the derivation, since the transformations are isometric, the kinetic energy term in the Hamil-tonianĤ remains invariant, andĤ commutes withD(S n ), the two operators consequently sharing the same eigenstates. We denote these eigenstates ψ aκ (r), with κ = (κ 1 , κ 2 , . . .). Hence we havê where λ(κ 1 ) is the eigenvalue ofD(S 1 ). Since electron density remains invariant under symmetry operations,

Operations
Examples of usage & notes S1(z) rotation bend tubes, wires, ribbons [4] S1(z) rotation and S2(z) translation bend membranes, slabs S1 joined rotation(z) + translation(z) twist nanotubes, wires, ribbons, DNA, simulate springs and coils [14-16] S1(x) and S2(y) two rotations around the same origin simulate spherical symmetry; solid and liquid membranes, such as mono-and multilayer graphene and lipid bilayers. (S1(x) and S2(y) commute approximately if rotation angles are small, and curvature can be treated as a local property.) S1(x), S2(y) rotations around different origins simulating arbitrary Gaussian curvature, like saddle structures (approximate treatment, like above) S1 translation(x)+reflection (yz) structures with repeating units . . . ABBA . . ., using AB unit; many waves can be simulated using half the wavelength[17] S1(x), S2(y) translations, S3(xy) reflection [+optional translation (x or y)] computational surface science; reflection doubles the surface slab thickness with half the number of atoms normal point group symmetries finite symmetric molecules and clusters By repeating the above steps for the remaining symmetry operations, we obtain a revised Bloch's theorem: in a potential, that is invariant in symmetry operations S n , the energy eigenstates ψ aκ at r and at r = S −n r differ by a phase factor exp(iκ · n), D(S n )ψ aκ (r) = ψ aκ (S −n r) = exp(iκ · n)ψ aκ (r). (5) This implies that wave functions only in one unit cellwhatever its shape-determine the electronic structure of the whole, extended system. You may recognize that the theorem above is nothing but Bloch's theorem, merely written with unconventional symbols. Indeed, many things remain as usual. Energy eigenstates ψ aκ automatically fulfill Eq.(5), once written in a revised version of Bloch basis, where ϕ µ (r) are local orbitals and n 1 = N is the number of unit cells. The Bloch basis gives Hamiltonian diagonal in κ, and similarly for overlap matrix elements. The total energy expressions remain the same, we only use a set of κ-points instead of k-points (extra symmetries reduce the set). [18] Because forces are calculated as parametric derivatives of the total energy, molecular dynamics works normally and energy is conserved; simulation cell dynamics, however, are different, and the concept of pressure needs redefinition. Finally, the theorem works with any electronic structure method, whether it is ab initio or not, whether it uses real-space grids or local orbitals (plane waves are tricky), or whether the approach is numerical or analytical.
Some things, however, do change in the revised Bloch's theorem. For bulk the periodic boundary condition is an approximation, whereas here some symmetries may form cyclic groups in reality, as in benzene. For cyclic groups the κ-point sampling is more restricted; in the above example of an achiral carbon nanotube the translational component κ 1 can be freely sampled between [−π, π] (because for S 1 periodicity is an approximation), but the rotational component accepts only the discrete values make the sampling of the components of κ coupled, for then we must have The connection between κand k-points is as follows. If S 1 is a translation along L 1 , then S 1 ψ aκ should give the same phase as T (L 1 )ψ ak , if ψ aκ and ψ ak are the same physical states. Hence exp(−iκ 1 ) = exp(−ik · L 1 ), and, in general, κ i = 3 j=1 L i j k j for i = 1, 2, 3. The formalism gives surprises, too. An atom can perform work on itself. This is because the total force on an atom exerted by its own periodic images-if rotations are involved-may differ from zero. Furthermore, a force F JI on atom I exerted by atom J may not be the counterforce to the force on atom J exerted by atom I, that is F JI = −F IJ ; Newton's third law appears invalid. These unorthodoxies are not bugs; remember that we simulate the whole extended system, and an atom I in the primitive unit cell is different from the atom I in a different unit cell-an artifact of atom indexing. Finally, note that if S i contains rotations, also local orbitals rotate; this is implicit in the operationD(S i )ϕ µ (r) in Eq. (6).
We implemented this formalism using local basis in the density-functional tight-binding software hotbit [19,20], and tested it with many finite and extended structures. We omit the details of the implementation here, and just comment on three things. First, the standard methods of electrostatics, like Ewald summation, are invalid since flexibility is required; we chose to use multipoles as they easily lend themselves for rotations and reflections. Second, the implementation can be done so that only the mappings are needed to build new symmetries; this requires just a couple of lines new code. Implementation generally is not hard, but it may be nontrivial for codes already build upon translational symmetry. Third, implementation has a negligible computational overhead compared to translational symmetry (see Table II). Certain manipulations take more time, but the most CPU-intensive parts remain as usual. So far our discussion has been abstract, but what can we do with the formalism in practice? While it may seem that we require a lot of symmetries, the main point of this Letter quite the opposite: we require less symmetries than before. Formalism enables simulating distorted materials, but also reduces computational costs for certain simulations. Selected examples of usage are shown in Table I. For example, one cost-reduction area is surface science, where less atoms are needed to simulate thick surface slabs. We believe more application areas can be discovered, once the new concepts are mastered. Now we leave the general discussion, and give one practical example of usage: we investigate twisted armchair graphene nanoribbons (AGNRs). [21] We choose this example for the existing literature, but also for the possibility to illustrate operations beyond standard chiral symmetry. Figure 1a shows a piece of an infinitely long 20-AGNR, with a twist χ = 360 deg/290.7Å= 1.2 deg/Å, within one conventional unit cell of 2720 atoms. The minimal unit cell, enabled by the new formalism, in turn, has 20 atoms (atoms A in Fig. 1b), and is accompanied by two symmetry operations: S 1 is L = 4.2Å translation, followed by L · χ = 5.04 degree rotation, and S 2 is L/2 = 2.1Å translation, followed by 182.52 degree rotation (= 180 deg+χ · L/2). The whole system can be built from one unit cell by S m1 1 S m2 2 with m 1 = 0, ±1, ±2, . . . and m 2 = 0, 1. Note that S 2 2 = S 1 and hence κ 2 = κ 1 /2 + πm 2 , while κ 1 is freely sampled. (While this was our choice for the symmetry operations, another, and equally sufficient, choice would have been to use only S 2 with m 2 = 0, ±1, ±2, . . ..) [22] Table II shows the wall-clock times for selected simulations. The simulations with χ = 0 show that the new formalism has no computational overhead compared to translation. The simulations with χ = 1.2 deg/Å show that finite twist affects simulation times with neither minimal nor chiral cells. The translational cell was too large for direct simulation, and the timing was estimated from the scaling law time∼(system size) 3 . Note that, by decreasing χ, the translational cell size-along with its timing-could easily be grown indefinitely. Point here is that twists even smaller than 1.2 deg/Å will be required to investigate the relevant physics of a 20-AGNR. Conventional quantum-mechanical simulation is practically impossible. We also remark that, given compatible κpoint samplings, energies and forces from different types of cells are the same within floating-point precision. Figure 1c shows AGNRs' energies as a function of twist. For ribbons wider than ∼ 12Å the energy is at minimum with non-zero χ = χ 0 ; ribbons twist spontaneously. This confirms earlier predictions by classical potentials (using thousands of atoms) [5,9] and finite el-  ement modeling [3]. The physical reason for twisting is the compressive edge stress that elongates edges with respect to ribbon's center. [5,9,23,24] The stress we get (1.7 eV/Å) agrees well with the stress (1.5 eV/Å) from previous density-functional calculations.  resides not only at the edge, but extends more into ribbon's center-a feature hard to reproduce by classical potentials. [9] This is also why we have no spontaneous twist for narrow ribbons. Ribbons ∼ 11Å wide have nearly zero torsion constant, and could be used in ultrasensitive torsion balances. Ultimately, very wide ribbons should show bifurcation into flat ribbons with ripples at the edges, but we won't discuss that here. [3,5] As argued in the abstract, the electronic properties ought to be investigated together with mechanical properties; this requires quantum mechanics. It is known that AGNRs have a gap due to the confinement of the finite width, and our gaps (inset in Fig. 2a) agree well with density-functional calculations of Ref. 21. But what happens to electronic structure when ribbons get twisted? Fig. 2a shows that twisting changes 10-AGNR's gap very little-this is generic for all AGNRs. The gap from πelectrons alone shows further that sp-rehybridization is negligible. Even the band structures of flat and twisted ribbons (Fig. 2b) are nearly identical. This suggests that, contrary to CNTs [25], the electronic properties of GNRs are remarkably robust against twisting.
To conclude, we hope to have illustrated how modest revision of Bloch's theorem enables versatile material distortions with quantum mechanics included, both numerically and analytically. However, excess emphasis on quantum mechanics causes undue discrimination of classical methods-the formalism works equally with classical force fields, finite element methods, or coarse-grained simulations, and equally when applied to, say, liquidphase cell membranes, fluid flow through bent pipes, or electron transport.
We acknowledge the Academy of Finland for funding, H. Häkkinen for support and the Finnish IT Center for Science (CSC) for computational resources.