Density functional theory method for twisted geometries with application to torsional deformations in group-IV nanotubes

We present a real-space formulation and implementation of Kohn-Sham Density Functional Theory suited to twisted geometries, and apply it to the study of torsional deformations of X (X = C, Si, Ge, Sn) nanotubes. Our formulation is based on higher order finite difference discretization in helical coordinates, uses ab intio pseudopotentials, and naturally incorporates rotational (cyclic) and screw operation (i.e., helical) symmetries. We discuss several aspects of the computational method, including the form of the governing equations, details of the numerical implementation, as well as its convergence, accuracy and efficiency properties. The technique presented here is particularly well suited to the first principles simulation of quasi-one-dimensional structures and their deformations, and many systems of interest can be investigated using small simulation cells containing just a few atoms. We apply the method to systematically study the properties of single-wall zigzag and armchair group-IV nanotubes, as they undergo twisting. For the range of deformations considered, the mechanical behavior of the tubes is found to be largely consistent with isotropic linear elasticity, with the torsional stiffness varying as the cube of the nanotube radius. Furthermore, for a given tube radius, this quantity is seen to be highest for carbon nanotubes and the lowest for those of tin, while nanotubes of silicon and germanium have intermediate values close to each other. We also describe different aspects of the variation in electronic properties of the nanotubes as they are twisted. In particular, we find that akin to the well known behavior of armchair carbon nanotubes, armchair nanotubes of silicon, germanium and tin also exhibit bandgaps that vary periodically with imposed rate of twist, and that the periodicity of the variation scales in an inverse quadratic manner with the tube radius.


Introduction
Over the past few decades, the synthesis and characterization of novel nanomaterials and nanostructures has blossomed into a major scientific and technological endeavor [1][2][3][4]. Such materials are usually associated with shapes and structures that are quite different from crystalline materials, and they often display properties that are radically distinct from the bulk phase. Consequently, a variety of computational techniques employing different physical theories have been developed over the years, to aid in their design and discovery [5][6][7][8][9].
A defining feature of the aforementioned class of materials is that they are of limited spatial extent along one or more dimensions. This often makes it possible to sustain unusual and/or large modes of deformation in such systems, without incurring material failure. Since a variety of material properties of nanostructures, including, e.g., optical, electronic and transport behavior are often strongly coupled to distortions in the material's structure, engineering the response of these systems through the application of mechanical strains constitutes an active and important area of scientific research today [10][11][12][13][14][15]. In particular, inhomogeneous strain fields -such as those associated with overall torsion (i.e., twisting) or flexure (i.e., bending) of the nanostructure, as well as those arising from localized deformations such as wrinkles or corrugations, have often been used to elicit fascinating electro-mechanical responses in such systems [16][17][18][19]. A persistent issue however, is that there appears to be a paucity of systematic and efficient computational techniques that can model these systems as they are undergoing such deformations, especially from first principles. We view the current contribution as an important step in addressing this gap in the literature and present a real-space formulation and implementation of Kohn-Sham Density Functional Theory (KS-DFT) that is suited to twisted geometries.
Systems associated with intrinsic twist are quite common among nanomaterials, with chiral carbon nanotubes [20], nanocoils [21] and inorganic nanoassemblies [22] constituting well known examples. Twisting is particularly relevant as a mode of deformation for quasione-dimensional systems such as nanotubes, nanoribbons, nanowires and nanorods [23], and can be an important route to engineering the properties of these materials through the imposition of strain. In particular, imposition of twist naturally gives rise to so-called helical potentials in achiral nanostructures, which can then cause these materials to display unusual transport properties and fascinating light-matter interactions [24]. Twisted geometries also have found relevance recently in the context of quasi-two-dimensional systems such as graphene bilayers [25][26][27][28], which are associated with strong electronic correlations and superconductivity, as well as the use of screw dislocations to engineer growth processes [29][30][31]. We anticipate that the simulation technique discussed in this work will have broad relevance to most of the materials systems described above, while being particularly consequential for the computational study of quasi-one-dimensional systems and their deformations, from first principles.
The vast majority of first principles calculations being carried out today use KS-DFT, as implemented using the pseudopotential plane-wave method [32][33][34][35]. While this is a powerful computational technique for the study of periodic systems (such as crystals) and their homogeneous deformations, it is fundamentally unsuitable for modeling systems subjected to inhomogeneous strain fields (such as those associated with bending or torsion), that break periodic symmetry. Indeed, modeling such systems by use of the plane-wave method can result in the use of uncontrolled approximations and/or performance and convergence (with respect to discretization parameters) issues that can render the calculations infeasible. For example, plane-wave calculations of a quasi-onedimensional system that is undergoing twisting (Figure 1) will usually involve making the system artificially periodic along the direction of the twist axis -thus resulting in a supercell containing a very large number of atoms, as well as the inclusion of a substantial amount of vacuum padding in the directions orthogonal to the twist axis, so as to minimize interactions between periodic images. Together, these conditions can make such calculations extremely challenging even on high performance computing platforms, if not altogether impractical.

Application of twist
Application of twist Axis of twist It has been pointed out in the literature however, that the aforementioned computational issues related to the study of twisted or bent nanostructures can be avoided by making use of the connections of such inhomogeneous strain states with non-periodic symmetries [23,[36][37][38][39][40][41][42][43][44]. Specifically, as long as edge effects are unimportant in a system under study, cyclic symmetries can be used to simulate bent nanostructures, while helical symmetries can be used to simulate systems with twist. A key ingredient for such an approach is the availability of efficient computational methods that can adequately handle such non-periodic symmetries. Following this line of thought, we have been developing systematic first principles simulation techniques suited to the study of systems with non-periodic symmetries [45]. In particular, we have developed ab initio methods that explicitly incorporate cyclic symmetries, and used this methodology to simulate bending in nanoribbons [46] and sheets of two-dimensional materials [47]. More recently, we have rigorously formulated and implemented a novel first principles computational technique that explicitly accounts for helical symmetries [48]. We view the present contribution as a follow up of this most recent development, and focus on the computational and application aspects of the simulation technique in this work, in contrast to our earlier contribution, which was largely concerned with the mathematical aspects. In particular, salient features of the current contribution are as follows. We present in this work a self-contained, intuitive derivation of the governing equations for systems associated with twisted geometries and make connections with helical symmetries, while also allowing for the possibility that such systems may have inherent cyclic symmetries. We describe the details of our computational strategy, including discretization choices in real and reciprocal space, numerical linear algebra issues and choice of eigensolvers. We touch upon specific aspects of our MATLAB based numerical implementation. We then discuss various features of the simulation method, including its convergence, accuracy, consistency, computational efficiency and parallel scaling properties. Finally, we apply the method to the study of torsional deformations of an important class of nanomaterials (i.e., nanotubes from Group IV of the periodic table 1 ) and investigate the electro-mechanical response of these systems. Notably, the present contribution subsumes our earlier work on KS-DFT for cylindrical geometries [47], and many of the results in that former contribution can be derived as special cases of the results presented here for twisted geometries (by considering simulations with zero twist). Together, the present contribution, and our earlier body of work extends symmetry adapted molecular dynamics and tight-binding based computational methods developed in the literature for studying bent and/or twisted nanomaterials, to the realm of first principles calculations.
The numerical technique described here employs finite difference discretization in helical coordinates 2 which allows us to set up a computational domain in an annular region of space. In turn, this enables us to carry out simulations of systems associated with twisted geometries, while employing small unit cells containing just a few atoms. With this setup in hand, we were able to carry out an extensive series of simulations involving zigzag and armchair nanotubes of carbon, silicon, germanium and tin, with radii approximately in the range of 1 to 3 nanometers. This enabled us to compare and contrast the properties of these different materials, and also allowed us to extend some well-known qualitative and quantitative features of the electro-mechanical properties of carbon nanotubes, to the broader class of Group IV nanotubes. We would like to point out that these studies would not have been possible without the use of a specialized computational method such as the one presented here. We anticipate that the rich repository of simulation data produced by our method can be utilized for the development of efficient, accurate, interpretable machine learning models [55], in the near future. 3 1 In modern IUPAC convention this group is also referred to as Group 14. Elsewhere, this group is also referred to as Group IVa or the Carbon group. 2 We are aware of chemistry literature based on Linear Combination of Atomic Orbitals (LCAO) methods [49][50][51][52][53][54], which have explored the use of helical and cyclic symmetries for studying nanostructures of interest. The connection of such symmetries with deformation modes in nanostructures does not appear to have been explored by these authors, as far as we can tell, and at any rate, these methods are quite distinct from the real space technique presented here. 3 After the submission of this manuscript, we were made aware of recent work [56] wherein the tech-

Formulation
In this section, we describe our formulation of Kohn-Sham density functional theory for twisted geometries. We first lay out the notation used in the rest of the paper. In what follows, e X , e Y , e Z will denote the standard orthonormal basis of R 3 . The Cartesian coordinates of a point p ∈ R 3 will be denoted as (x p , y p , z p ), i.e., x = x p e X +y p e Y +z p e Z . The corresponding helical coordinates (introduced later in Section 3.1) and cylindrical coordinates of the point will be denoted as (r p , θ 1 p , θ 2 p ) and (r p , ϑ p , z p ) respectively. The coordinates of a generic point will be denoted as (x, y, z), (r, θ 1 , θ 2 ) and (r, ϑ, z) in Cartesian, helical and cylindrical coordinates respectively. Vectors and matrices will be denoted in boldface, with vectors typically denoted using lower case letters (e.g., p) and matrices using uppercase (e.g. Q). The symbol · will be often used as a generic placeholder instead of specifying a variable explicitly (e.g. f (·) instead of f (x) or f (y)). The notation L 2 (Ω) will be used to denote the space of square integrable functions over a domain Ω. The inner product over such a space will be denoted as ·, · L 2 (Ω) . An overbar will be used to denote complex conjugation (e.g. f (x)). Finally, |·| will be used to denote the absolute value of a scalar, and · will be used to denote the norm of a vector or function.

System specification: Computational domain, atomic configuration and symmetries
We consider a nanostructure aligned along e Z , the axis of twist, as the prototypical system of interest. In order to avoid quantum finite-size effects and/or mechanical constraints at the edges due to the imposition of twist [40,57] , we will assume that the structure is infinite in extent along e Z . For the sake of simplicity, we will also assume that the structure is of limited spatial extent along e X and e Y , i.e., it is a quasi-onedimensional system. The large majority of nanomaterials for which twisted geometries might be relevant as deformation modes, are included within the scope of the above set of assumptions. These conditions imply that the system can be embedded in a cylinder with axis e Z (or annular cylinder, if the system is tubular), of infinite height and finite radius, and we will refer to this region of space as the global simulation domain.
For most quasi-one-dimensional systems of interest, the infinite extent along e Z is related to periodicity along this axis. Additionally, for many such systems, including for example, the tubular structures considered in this work, there may be rotational symmetries about the same axis. Let the atoms of the untwisted structure have positions: The above assumptions on periodicity and rotational symmetry imply that there is a periodic group consisting of translations along e Z : niques presented by us here as well as our earlier contribution [48] have been implemented into an efficient C/C++ framework. a cyclic group of order N about e Z (consisting of rotations through multiples of the angle and a finite collection of points: such that the entire structure S untwisted can be described as the action of the composite group: on the points in P, i.e., In the above equations, a symbol of the form Q | t denotes an isometry with rotation Q ∈ SO(3) and translation t ∈ R 3 . Its action on a point x ∈ R 3 can be written as: Additionally, R nΘ denotes the following rotation matrix with axis e Z : I denotes the identity matrix and 0 denotes the zero vector. The scalar 0 < τ < ∞ is the fundamental period of the group G periodic . We will refer to the points in P as the simulated atoms. We will use Z k to denote the valence nuclear charge of the simulated atom located at position r k . Now let us suppose that the structure S untwisted is subjected to a uniform twist of 2πα radians per τ bohr along the axis e z , so as to result in the structure S twisted with the atomic positions: Since we are dealing with structures that extend to infinity along e Z , we may obtain the deformed (twisted) configuration by prescribing a mapping of the form q = R2παzp τ p, to the undeformed one. Here, α ∈ [0, 1) is a scalar twist parameter, τ can be re-identified as the pitch of the twist, and β = 2πα τ , is the rate of twist. Furthermore, denotes a rotation matrix with axis e Z for which the (twist) angle depends on the coordinate along e Z . At the atomic level, this implies [23,36,48] that the deformed structure may be obtained from the undeformed one by replacing the group of translations G periodic used to generate S untwisted , by a group of screw transformations (or helical isometries), i.e.: Here R 2πmα denotes the following rotation matrix with axis e Z : In other words, by replacing the composite group G untwisted with: we may generate the structure with the prescribed amount of twist as: In the above equations, R (2πmα+nΘ) denotes the following rotation matrix with axis e Z : Note that in this formulation, the structure continues to maintain its cyclic symmetries even after twisting. Also note that the formula in eq. 14 (and similarly, eq. 6) is meant to be species preserving in the sense that an atom in the simulated set P has the same atomic number as its images under the isometries in G untwisted (or G untwisted ). 4 Also note that by virtue of the above definitions, the group G twisted serves as a physical symmetry group for the structure S twisted in the sense that the action of any Υ ∈ G twisted on all the points in S twisted leaves it invariant (and similarly for G untwisted and S untwisted ). The group G twisted will play a central role in the rest of this work. Note that this group subsumes the group G untwisted in the sense that the latter can be recovered by simply setting α = 0 in the former. In what follows, we will simplify notation a bit and simply use G to denote this group. Further, we will use the notation: to denote group elements from G. The action of Υ m,n on a generic point in space is to rotate it about axis e Z by angle 2πmα + nΘ while also translating it by mτ along the same axis.
In subsequent sections, we will describe how the Kohn-Sham problem for the entire twisted structure as posed on the global simulation domain, can be appropriately reformulated as a problem over a fundamental domain (or symmetry adapted unit cell ), such that only the simulated atoms and the symmetry group G are involved in the resulting equations. This symmetry adapted computational domain has to be a regular region of space with sufficiently smooth boundaries that encompasses the simulated atoms and can be used to tile the global simulation domain by the action of the group G. Furthermore, this region should be minimal in the sense that the above tiling operation should not produce intersecting volumes. In the context of the twisted tubular structures considered in this work, if the simulated atoms have radial coordinates lying between R in and R out , the following region serves as an appropriate fundamental domain (expressed using cylindrical coordinates): The boundaries of the fundamental domain defined above can be expressed as: Here ∂R in and ∂R out denote boundaries related to the radial direction (i.e., the surfaces r = R in and r = R out respectively), ∂ϑ 0 and ∂ϑ Θ denote (z-dependent) bounding surfaces related to the angular direction (i.e., ϑ = 2παz τ and ϑ = 2παz τ + Θ respectively), and finally, ∂Z 0 and ∂Z τ denote boundaries related to the e Z direction (i.e., the surfaces z = 0 and z = τ respectively). Note that for no applied twist, the region D is simply an annular cylindrical sector, i.e., and the boundaries ∂ϑ 0 and ∂ϑ Θ are then vertical surfaces perpendicular to the e Y − e Z plane. Figure 2 shows two views of the fundamental domain used for the simulations used in this work and also highlights the boundaries described above. In what follows, we will formulate suitable versions of the equations of Kohn-Sham theory as posed on the simulation cell D and also elaborate on the conditions that have to be applied on the bounding surfaces that make up ∂D. Our derivation of the governing equations presented here is largely heuristic, and a more nuanced, mathematically rigorous discussion is available in [48].
2.2. Governing equations 2.2.1. Helical Bloch theorem and block-diagonalization of Hamiltonian As described above (eq. 14), the atomic positions of the twisted structure can be described as the orbit of a discrete group of isometries (i.e., the group G). Due to the presence of such symmetries in the system, it follows under fairly general hypotheses [45][46][47][48] that the ground state electron density for such a system is invariant under the same symmetry group. Furthermore, the Kohn-Sham Hamiltonian for the system commutes with the symmetry operations of the group [58,59]. Consequently, the eigenstates of the Hamiltonian can be labeled using irreducible representations of the group G, and they transform under action of the group in the same manner as the irreducible representations themselves do [45,48,58,59]. Since the group G is Abelian, results from group representation theory [60,61] imply that the complex irreducible representations are one dimensional. These are the so called complex characters of G, which, keeping in mind that G is the direct product of the groups G helical and G cyclic , can be expressed as (for m ∈ Z, n ∈ {0, 1, 2, . . . , N − 1}): In other words, for each value of η ∈ and ν as defined above, the characterζ ∈ G is a complex valued map on the group 5 , that assigns the value e 2πi mη+ nν N to the group element Υ m,n ∈ G. Since any characterζ ∈ G can be labeled using the pair (η, ν), these can be also used to label the eigenstates of the Kohn-Sham Hamiltonian, and other quantities related to its spectrum. Accordingly, we will use λ j (η, ν), ψ j (x; η, ν) and g j (η, ν) to explicitly indicate the labels for the eigenvalues, the eigenvectors, and the electronic occupations, respectively. Collections of the eigenvalues, eigenvectors and electronic occupations will be denoted using Λ, Ψ and G respectively, i.e.: ; ν ∈ 0, 1, 2, . . . , N − 1 ; j = 1, 2, . . . , ∞ , ; ν ∈ 0, 1, 2, . . . , N − 1 ; j = 1, 2, . . . , ∞ , Mathematical properties of the characters and the above discussion lead to a number of important considerations that are worth mentioning at this point. First, as a consequence of the orthogonality relations obeyed by the characters [60,61] the eigenstates associated with distinct characters are orthogonal to each other. This can be used to cast the Hamiltonian (which commutes with the symmetry operations in G) in a symmetry adapted basis [59], such that it appears block-diagonal [45,48]. Since the blocks associated with distinct characters can be dealt with independently of each other and are of reduced dimension compared to the full Hamiltonian (within any finite dimensional approximation, e.g.), this implies that the problem of diagonalizing the Hamiltonian is greatly simplified. Second, the fact that the eigenstates of the Hamiltonian transform under symmetry operations in the same manner as the characters, implies that they obey a Helical Bloch theorem [45,48,57,62], i.e., for any Υ m,n ∈ G: or equivalently: These relations can be used to deduce the conditions that need to be applied to the boundary surfaces of the fundamental domain while formulating the Kohn-Sham problem. Finally, in order to write down quantities that depend on all eigenstates cumulatively, we need to account for contributions from each ζ ∈ G. This amounts to integrating the eigenstate dependent quantities against a suitable integration measure over G, i.e., by forming sums of the form , along with integrals in η. As an example, if we intend to compute the sum of the occupation numbers over all the electronic states in the system, we need to evaluate: Here and henceforth, I is used to denote the set − 1 2 , 1 2 .

Electronic free energy functional and Kohn-Sham equations for twisted structure
In what follows, we will consider the (twisted) system of interest to be one in which the effects of spin can be ignored, and for which the electronic temperature is set at T e . This implies that the electronic occupations can be expressed in terms of the Kohn-Sham eigenvalues as: with f Te ·) denoting the Fermi-Dirac function, i.e., Here λ F and k B denote the system's Fermi level and the Boltzmann constant respectively. In order to motivate the correct form of the various terms of the governing equations for the twisted structure, we will often refer to the simpler, more well known expressions of these quantities for finite (or isolated) systems. We will denote these finite system relevant quantities (scalar fields, energies, etc.) with a • superscript. For a finite system [47,63], the electron density can be expressed in terms of the Kohn-Sham eigenvectors and the electronic occupations as: Following the discussion above, this expression has to be modified for our case as: Note that the factor of 2 in the expressions above is due to ignoring electronic spin. Further note that due to the Helical Bloch conditions obeyed by the Kohn-Sham eigenvectors (eq. 23), the expression above is invariant under the symmetry operations in G, as is required of the ground state electron density.
• Electronic free energy: To derive the governing equations of Kohn-Sham theory for our system, we take recourse to an energy minimization approach [47,48,64]. The relevant quantity in this case, since the system is of an extended nature, is the ground state electronic free energy per unit fundamental domain. We denote this quantity here as F(G, Ψ, P, D, G) to emphasize its dependence on the electronic occupation numbers, the eigenstates, the positions of the simulated atoms, the fundamental domain and the symmetry group G. Within the pseudopotential [8,65] and Local Density Approximations [66], we may express it as: The terms on the right-hand side of the above expression represent (per unit fundamental domain) the kinetic energy of the electrons, the exchange correlation energy, the nonlocal pseudopotential energy, the electrostatic energy and the electronic entropy contribution, respectively. We now elaborate on each of these quantities.
• Kinetic energy: The first term on the right hand side of the expression above is the electronic kinetic energy per unit fundamental domain. For an isolated system (placed in R 3 ), this term can be expressed [47,63] in terms of the Kohn-Sham eigenstates and the occupations as: For the system at hand, this is modified to read: • Exchange-correlation energy: The second term represents the exchange correlation energy per unit fundamental domain and is expressible using the Local Density Approximation (LDA) [66] as: Note that the above formulation does not preclude the use of more sophisticated exchange correlation functionals such as the Generalized Gradient Approximation [67]. Since the use of such functionals has little bearing on the subsequent discussion, we do not consider them further in this work.
• Nonlocal pseudopotential energy: The third term on the right hand side of eq. 29 represents the nonlocal pseudopotential energy per unit fundamental domain. For a finite system consisting of M • atoms located at the points {r • k ∈ R 3 } M • k=1 , the non-local pseudopotential operator in Kleinman-Bylander form [68] can be written as: in terms of the projection functions {χ k,p (·; r k )} N k p=1 and the corresponding normalization constants {γ k,p } N k p=1 associated with the k th atom (located at y k ). The nonlocal pseudopotential energy in that case has the form: To obtain the analogous expression for the twisted structure, we consider the contributions from the atoms located within the fundamental domain and all the electronic states in the system [48] to get the nonlocal pseudopotential energy per unit fundamental domain as: Here, the overlaps of the orbitals with the atom centered projectors are carried out over the global simulation domain C, since the latter can have support extending beyond the fundamental domain. With the aid of the Helical Bloch Theorem (eq. 23) and by using the properties of the projection functions χ k,p , the integral implicit in the above expression can be reduced to the fundamental domain [47,48], so that a more computationally convenient expression for the nonlocal pseudopotential energy per unit fundamental domain reads as: The functionsχ k,p in the equation above can be expressed as: • Electrostatic interaction energy: The fourth term on the right hand side of eq. 29 represents the electrostatic interaction energy per unit fundamental domain. This includes the Coulombic attraction between the electrons and the nuclei, as well as the mutual repulsion between the electrons themselves. To express this term, it is useful to introduce the net electrostatic potential Φ, which also appears in the Kohn-Sham equations (as part of the effective potential). To see how this can be done, we consider first a finite system placed in R 3 , with nuclei located at the points {r • k ∈ R 3 } M • k=1 . For this example system, the net electrostatic potential Φ • , can be expressed in terms of the total charge of the (finite) system as: Here, ρ • represents the electron density and b finite represents the total nuclear pseudocharge. The latter can be expressed in terms of the individual nuclear pseudocharges Note that for each atom, the pseudocharge (typically a smooth, radially symmetric, compactly supported function) integrates to its valence nuclear charge, i.e., The connection between the potential Φ • and the electrostatic interaction energy is that we may express the latter as: and the scalar field Φ • which attains the maximum in the above problem is precisely the one presecribed using eq. 38. Note that the constant term E • sc (r • 1 , r • 2 , . . . , r • k ) is added as a correction for self-interactions and possible overlaps of the nuclear pseudocharges [69]. It is independent of Φ • and does not play a role in the above optimization problem.
With the above discussion in mind, we may now introduce the net electrostatic potential for the twisted structure using the electron density (eq. 28) and the net nuclear pseudocharge associated with the system, in a manner that is analogous to eq. 38, i.e., The net nuclear pseudocharge at any point in the global simulation domain can be expressed using the pseudocharges of the atoms in the fundamental domain as: Note that since the electron density is group invariant, as is the net nuclear pseudocharge (by construction), the total electrostatic potential for the twisted structure is group invariant as well. Thus, it suffices to compute this quantity over the fundamental domain, in addition to specifying boundary conditions that are consistent with the group invariance of the function. Following eq. 41, we now write the electrostatic interaction energy per unit fundamental domain as: The scalar field Φ which attains the maximum in the above problem, is the same one specified in eq. 42. The constant (i.e., Φ-independent) term E sc (P, G, D) accounts for self-interaction corrections and possible overlaps between pseudocharges. We omit the details of this term here for the sake of brevity, and cite references [46,63,69] for relevant details.
• Electronic entropy: Finally, the last term on the right hand side of eq. 29 deals with the contribution of the electronic entropy to the free energy. Using Fermi-Dirac smearing, for a finite system at electronic temperature T e , the electronic entropy can be represented : Analogously, the corresponding term for the twisted structure reads as: • Kohn-Sham Equations: With the expressions for the various energy terms in place, we write the electronic ground-state energy for the twisted structure as the following minimization problem: with the added constraints that: 1. the orbitals in Ψ are helical Bloch states, namely, they obey eq. 23 and are orthonormal over the fundamental domain for each ζ ∈ G, i.e.: and, 2. the number of electrons per unit fundamental domain is a fixed number, i.e., The Euler-Lagrange equations corresponding to the above variational problem are the Kohn-Sham equations for the twisted structure, as posed on the fundamental domain. For j ∈ N, η ∈ I and ν = 0, 1, . . . N − 1, we may express them as: with H KS denoting the Kohn-Sham operator, i.e.: Here, V xc denotes the exchange correlation potential: Φ (as introduced in eq. 42) denotes the net electrostatic potential arising from the electrons and the nuclear pseudocharges, and obeys the Poisson equation: while V nl denotes the non-local pseudoptential operator (specifically, its (η, ν) component), and can be expressed in terms of the functionsχ k,p (introduced in eq. 37) as: Note that the use of eq. 53 in lieu of eq. 42 is preferable for practical calculations since computationally inconvenient non-local integrals that appear in the latter equation are avoided [47,[69][70][71]. Together, eqs. 50 -54, along with eq. 48 and 49 form the governing equations for our system and need to be solved self-consistently.

Boundary Conditions
The unknown fields in the governing equations above are the orbitals ψ j (·; η, ν) ∈ Ψ and the electrostatic potential Φ. Since these fields are posed on the fundamental domain D, we need to augment the governing equations with boundary conditions on the surfaces that make up ∂D. By using the conditions in eq. 23 on the orbitals, and observing that the symmetry operation Υ 1,0 = R 2πα | τ e Z maps ∂Z 0 to ∂Z τ , while the operation Υ 0,1 = R Θ | 0 maps ∂ϑ 0 to ∂ϑ Θ , we arrive at: Concurrently, since the net electrostatic potential is invariant under all symmetry operations in G, it obeys the boundary conditions: The above equations leave the boundary conditions on the surfaces ∂R in and ∂R out unspecified. As far as the wavefunctions are concerned, we may enforce Dirichlet boundary conditions on these surfaces, by appealing to the decay of the electron density along the radial direction [47,48]. This gives us: On the other hand, the electrostatic potential Φ may not decay to zero quickly along the radial direction. Therefore, it is more prudent to set Φ(x ∈ ∂R in ) and Φ(x ∈ ∂R out ) by direct evaluation of eq. 42 by using a modified version of the Ewald summation technique [72]. In practical calculations however, this correction may be sometimes ignored [48] in favor of Dirichlet boundary conditions on those surfaces.

Other quantities of interest at self-consistency
At the end of the self consistent field iterations, a number of other quantities may be computed from the converged electronic states. For instance, we may obtain a more accurate estimate (i.e., one that is less sensitive to self-consistency errors) of the Kohn-Sham ground state electronic free energy (per unit fundamental domain) by using the Harris-Foulkes functional [73,74] instead of eq. 29. This can be written for the twisted structure, using quantities expressed over the fundamental domain as: Note that the first term on the right hand side of the above equation is the electronic band energy.
For ab initio molecular dynamics or structural relaxation calculations, atomic forces need to be calculated. The Hellmann-Feynman forces on the atom located at r k in the fundamental domain can be computed about the ground-state as: Note that since the forces are derivatives of a free energy which is invariant with respect to the symmetry operations in G, it follows that the force on an atom Υ m,n •r k located outside the fundamental domain can be evaluated in terms of the force on its counterpart in the fundamental domain as (R 2πmα+nΘ ) −1 f k [36]. Thus, to perform structural relaxations on the twisted structure, it suffices to concentrate on the atoms in the fundamental domain and drive their forces to zero. Finally, the electronic density of states which often offers useful information about the electronic properties of a material under study, can be computed at an electronic temperature T e as [75]: with f Te (·) denoting the derivative of the Fermi-Dirac function.

Implementation
We now discuss different numerical and computational aspects of the implementation of the above methodology.

Use of helical coordinates
The equations in Section 2 above are expressed in a manner that do not make any explicit reference to a coordinate system. For numerical implementation purposes however, it is useful to introduce a coordinate system that is commensurate with the geometry of the twisted structure and the symmetries of the system. The helical coordinate system, introduced in [45,48] is well suited for these purposes. If a point p in the global simulation domain C has Cartesian coordinates (x p , y p , z p ) and cylindrical coordinates (r p , ϑ p , z p ), then the corresponding helical coordinates (r p , θ 1 p , θ 2 p ) are defined as: 17 The helical coordinates reduce to the usual cylindrical coordinates when the twist parameter of the system is 0 and the pitch τ is set to unity. The inverse relations: map the helical coordinates of p to their Cartesian counterparts. The coordinate transformations introduced above can be used to map the curvilinear coordinate system associated with the twisted structure, to a rectilinear one in which computations are simpler to set up. Specifically, the relations in eq. 64 above map the cuboid (R in , R out ) × (0, 1) × (0, 1/N) to the fundamental domain D. In particular, the bounding surfaces of the fundamental domain can be described in helical coordinates as r = R in (for ∂R in ), r = R out (for ∂R out ), θ 1 = 0 (for ∂Z 0 ), θ 1 = 1 (for ∂Z τ ), θ 2 = 0 (for ∂ϑ 0 ) and θ 2 = 1/N (for ∂ϑ Θ ). Furthermore, the symmetry operation Υ m,n maps the helical coordinates of a point p from (r p , θ 1 p , θ 2 p ) to (r p , θ 1 p + m, θ 2 p + n N ). In order to express the equations in Section 2.2 in helical coordinates, we need the the Laplacian operator, the Cartesian gradient and the integral of a function (over the fundamental domain) expressed in helical coordinates. For a function f (r, θ 1 , θ 2 ) these take the form [48]: Upon expressing the Kohn-Sham orbitals as ψ j (r, θ 1 , θ 2 ; η, ν), the above expressions allow the governing equations and boundary conditions to be expressed in helical coordinates exclusively. For numerical implementation purposes however, it is convenient to work with functions that are completely invariant under symmetry operations instead of being invariant upto a Bloch phase, as the orbitals are. To this end, we write: where the functions u j (r, θ 1 , θ 2 ; η, ν) are group invariant. In terms of these auxiliary functions, the governing equations over the fundamental domain are: ρ(r, θ 1 , θ 2 ) = 2 The boundary conditions 6 are:

Approximation of infinite series in governing equations
The governing equations as posed above, contain series sums over infinite numbers of terms which need to be truncated for the purposes of numerical implementation. Such infinite sums not only appear explicitly while summing over an infinite number of electronic states (eqs. 71, 72), but also implicitly in the calculation of quantities such as the net pseudocharge (eqs. 70,43) and the nonlocal pseudopotential operator (eqs. 69,54,37). We now describe our strategies for dealing with such quantities.
In order to truncate sums involving an infinite number of electronic states, we may assume -as is commonly done in the literature [77,78], that the electronic occupation numbers decay to zero beyond the lowest N s electronic states. In effect, this implies that sums over the index j in equations 71 -72 run from 1 to N s (instead of 1 to ∞), and that a set of N s eigenvalue problems for each value of η and ν, as posed in eq. 69, have to be considered. In practical calculations when the electronic temperature is less than a few thousand Kelvin, the number of states N s can be related to the number of electrons per unit cell N e as N s = c s × N e 2 , with the constant c s chosen to be between 1.05 and 1.20 [78]. The infinite sums involved in calculation of the net pseudocharge and the non-local pseudopotential operator both arise due to summations over individual atoms in the fundamental domain, as well as their repeated images under the group G (eqs. 43,54,37). However, we observe that the functions being summed in these cases are always centered about the atoms in question, and they have the property of being compactly supported in a small spherical region of space around the atom (i.e., the functions b k (·) in eq. 43 and χ k;p (·) in eq. 37 all have this property). Thus, the contribution of such terms to the fundamental domain is zero beyond a few terms of the series expressed in eqs. 43 and 37, and the infinite summations in these expressions can be reduced to a set of values of m and n that are near zero. 7

Discretization Strategy
The equations above need to be discretized in real space (i.e., over the fundamental domain D) as well as in reciprocal space (i.e., over the set G). We now describe our strategies for addressing each of these issues.

Real space discretization of the fundamental domain
We use a higher order finite difference scheme [47,48,63,64,[79][80][81][82] for real space discretization. Since helical coordinates have the property of "unwrapping" the fundamental domain D to a cuboid, a convenient meshing of the computational domain can be attained by choosing equispaced points in the r, θ 1 and θ 2 directions. Accordingly, we choose N r , N θ 1 and N θ 2 grid points along these directions (respectively), and observe that the corresponding mesh spacings h r , h θ 1 , h θ 2 satisfy: We will often refer to the quantity h = Max. h r , τ h θ 1 , 2π R in +Rout 2 h θ 2 as the overall mesh spacing for a particular level of discretization. We index each finite difference node using a triplet of natural numbers: and we use f (i,j,k) to denote the value a function f at the grid point i, j, k. The grid point with indices (i, j, k) is associated with the radial coordinate To discretize equations 69 -72 using the finite difference scheme, we require expressions for first and second order derivatives, as well as a quadrature rule to compute integrals over the fundamental domain (e.g., to evaluate the action of V nl on a given function). The expressions for the first order derivatives are: The expressions for the second order derivatives are: , , In the above expressions, n o denotes half the finite difference order, s denotes r, θ 1 or θ 2 , and the finite difference weights w second p,s and w first p,s can be expressed as [83]: Finally, the expression for evaluating integrals over the fundamental domain is:

Reciprocal space discretization
As is evident from the governing equations, many quantities of interest (including the electron density, for example) involve accumulating sums from each of the characters ζ ∈ G. Since this is equivalent to computing sums of the form 1 N N−1 ν=0 and integrals over the set I, we need a suitable scheme for discretizing the set G. Accordingly, we perform quadratures over the set G using: In the above expression, in accordance with the Monkhorst-Pack scheme [84], the quadrature nodes η b are equi-spaced, while the corresponding quadrature weights w b are uniform. Effectively, the above scheme discretizes the set G using N K = N η × N representative reciprocal space points. By use of time reversal symmetry [47,48,85], it is possible to reduce the number N K by a factor of 2, which helps in cutting down computational wall time in practical calculations.

Solution strategies for the discretized equations and MATLAB implementation
The governing equations for the twisted structure represent a set of coupled nonlinear eigenvalue problems. We use self consistent field (SCF) iterations accelerated via Periodic-Pulay extrapolation [86] to solve them in this work. The total effective potential (i.e., V xc + Φ) is used as the mixing variable. Solution of the Poisson equation associated with the electrostatic field (eq. 70) is carried out using the Generalized Minimal Residual method (GMRES) [87], and an incomplete LU factorization based preconditioner [88] is used to accelerate convergence of the GMRES iterations. Solution to eq. 72 is carried out using a nonlinear equation root finder [89].
As a consequence of the discretization choices and other simplifications outlined previously, there are N K linear eigenvalue problems, each of dimension N D , that have to be solved on each SCF iteration step. Furthermore, for each of these eigenvalue problems, the lowest N s eigenstates have to be determined via a suitable diagonalization process. Due to our use of finite differences, the discretized Hamiltonian operators (at each value of η and ν) are non-Hermitian, even though the original infinite dimensional operators they represent are not. This is a well known issue that arises while approximating differential operators such as the Laplacian in curvilinear coordinates using finite differences [46,47,90]. In practice, this issue is mitigated by a combination of factors. First, as the mesh spacing h is made finer, and/or the degree of the finite difference discretization n o is made larger, the discretized operators approximate their infinite dimensional counterparts more closely. Consequently, the discretized operators become more Hermitian (i.e., the norm of the difference between the operator and its Hermitian conjugate goes to zero), the eigenvalues have small imaginary components, and conventional iterative methods for obtaining the spectrum of a sparse symmetric Hamiltonian [91][92][93] can be employed for diagonalization. Indeed, for the discretization parameters used to produce the results in this work, the imaginary parts of the Hamiltonian eigenvalues are small enough that they can be ignored without any deleterious effects on the stability or quality of the simulations. Second, by choosing eigensolvers that can handle non-Hermitian problems in a robust manner, even calculations involving relatively coarse meshes (i.e., for which the Hamiltonian is well conditioned, but might have some eigenvalues with nonvanishing imaginary parts), or problems with poorly conditioned Hamiltonian matrices (which can arise if a system with a large amount of prescribed twist is being studied) can be performed.
Keeping the above factors in mind, our implementation employs a combination of the Generalized Preconditioned Locally Harmonic Residual (GPLHR) method [94], as well as iterative diagonalization based on Chebyshev polynomial filtered subspace iterations (CheFSI) [91,95,96]. Due to the ability of GPLHR to employ preconditioners (based on incomplete LU factorization, e.g.), the method can be particularly effective in handling poorly conditioned Hamiltonian matrices -i.e., for problems in which the CheFSI method tends to use relatively large polynomial filter orders. For such problems, we have also observed that GPLHR generally requires fewer iterations to reach SCF convergence, when compared to CheFSI, and between 5 to 8 iterations of the method are sufficient per SCF step. Nevertheless, for the systems considered in this work, we found that Chebyshev polynomial filter orders in the range 60 to 80 were adequate in guaranteeing stable, well converged simulations, and in this scenario the CheFSI method generally required shorter wall-times-to-solution overall. Thus, for the bulk of the simulations presented in this work, CheFSI was the method of choice. We show examples of the SCF convergence behavior for two example systems using CheFSI and GPLHR in Figure 3.

SCF residual
Si armchair nanotube, CheFSI Si armchair nanotube, GPLHR C zigzag nantube, CheFSI C zigzag nantube, GPLHR Figure 3: Examples of SCF convergence using CheFSI and GPLHR methods in Helical DFT. The armchair silicon nanotube has radius 0.96 nm, and was subjected to a twist of 5.67 degrees/nm. The zigzag carbon nanotube has radius 0.70 nm, and was subjected to a twist of 4.27 degrees/nm.
We have implemented the above methods and algorithms in a computational package called Helical DFT. The current version of the code is largely written in MATLAB [97], with certain key routines (including Hamiltonian matrix-vector products, sections containing multiple nested loops, etc.) written in C to alleviate speed and/or memory footprint issues. The code makes use of MATLAB's vectorization capabilities, and achieves parallelization by performing diagonalization of the Hamiltonian for different values of η and ν simultaneously over multiple computational cores. Helical DFT is capable of performing structural relaxation by use of the Fast Intertial Relaxation Engine (FIRE) algorithm [98] as well as ab initio molecular dynamics simulations by use of a velocity Verlet integrator [99].

Computational Platform
All simulations involving Helical DFT were run using a dedicated desktop workstation (Dell Precision 7920 Tower) or single nodes of the Hoffman2 cluster at UCLA's Institute for Digital Research and Education (IDRE). The desktop workstation has an 18-core Intel Xeon Gold 5220 processor (24.75 MB cache, running at 2.2 GHz), 256 GB of RAM and a 1 TB SATA Class 20 Solid State Drive (SSD). Each compute node of the Hoffman2 cluster has two 18-core Intel Xeon Gold 6140 processors (with 24.75 MB cache, running at 2.3 GHz), 192 GB of RAM and local SSD storage. MATLAB version 9.7.0 (R2019b) was used for the simulations. Compilation of C language routines was carried out using MinGW (on the workstation) and GCC (on the Hoffman2 nodes) software suites. Interfacing between MATLAB and C language routines was carried out by means of MATLAB's MEX and Coder frameworks, while parallelization was achieved by use of using MATLAB's Parallel Computing Toolbox.

Simulation Parameters
We used an SCF iteration convergence tolerance of 10 −6 in the total effective potential (relative residual). The Periodic Pulay mixing scheme [86] used a history of 7 iterations, the mixing parameter was set at 0.2, and Pulay extrapolation was performed on every alternate SCF step. GMRES iterations for the Poisson problem was carried out till the residual dropped below 10 −9 on every SCF step. We employed an electronic temperature of T e = 315.77 Kelvin in the Fermi-Dirac function (this corresponds to about 1 milli-Hartree of smearing), and included 2 extra states at each value of η and ν to accommodate fractional occupancies. We used Troullier-Martins norm conserving pseudopotentials [65] and Perdew-Wang parametrization [100] of the Local Density Approximation [66]. We used a 12 th order finite difference discretization scheme (i.e., n o = 6 in eqs. 77,78,79) and diagonalization via CheFSI used filters of order 60 to 80. Determination of spectral bounds within the CheFSI method used MATLAB's eigs function [101] with a relatively loose tolerance of 10 −2 . For the nanotube simulations described here, we ensured a gap of 10 to 11 Bohrs between the atoms located within the fundamental domain, and the boundary surfaces ∂R in and ∂R out , in order for the electron density and the wavefunctions to decay sufficiently in the radial direction 8 . Real space and reciprocal space discretization parameters were chosen on a case by case basis, as described later.

Materials Systems: Group IV Nanotubes
Nanotubes and other similar systems are particularly well suited for study using the methods described in this work. We choose single walled nanotubes of carbon, silicon, germanium, and tin as materials systems for investigation here. These systems are used for carrying out numerical validation studies, and due to their technological importance, also for gaining insights into their properties by the use of our method. Such one-dimensional nanostructures from Group IV of the periodic table can be described in terms of a "roll-up" procedure [102], starting from their two-dimensional sheet counterparts (i.e., graphene, silicene, germanene and stanene). We collectively refer to these oneand two-dimensional materials as X (X = C, Si, Ge, Sn) nanotubes, and Xenes, respectively. Both these classes of materials have been intensely studied in recent years through both experimental and computational methods, due to their association with fascinating materials properties [103-107, 107-109, 109-142]. In particular, the electronic properties of deformed carbon nanotubes have received extensive treatment in the literature through theoretical and computational means [62,135,[143][144][145][146][147][148][149][150]. Although a few computational studies on the electronic structure of the larger class of Group IV nanotubes are also available [47,132,[151][152][153], as far as we can tell, this work is the first to investigate from first principles, the behavior of these materials under torsional deformations, and to extend some well known results for carbon nanotubes to the broader class of Group IV nanotubes.
By using the roll up construction on the Xene sheets (see Figure 4), we can represent untwisted tubes using just four atoms in the fundamental domain [23,36,47], and a twist can be prescribed on the system by choosing a non-zero value of α. Depending on the direction of rolling, the untwisted tubes can be classified as armchair or zigzag, and the fundamental period τ of the untwisted tubes in these cases are √ 3 a and 3 a , respectively, with a denoting the (planar) interatomic distance among the X atoms. Furthermore, the cyclic group order N can be expressed in terms of the nanotube radius via the relation N L = 2πR avg. . Here L = 3 a and √ 3 a, for armchair and zigzag cases, respectively, and R avg. denotes the average radial coordinate of the atoms in the fundamental domain. For subsequent simulations, we adopted the values of the parameter a, as well as the out of plane buckling parameter δ, as reported in [47]. We include the values of the parameters in Table 1

Convergence, accuracy and efficiency studies
We begin with a discussion of the convergence properties of our numerical implementation with respect to discretization parameters. We choose armchair nanotubes of carbon (radius = 1.07 nm, N = 16), silicon (radius = 0.97 nm, N = 9), germanium 9 To compute these parameters, the relaxed ground state structure of the Xene sheets (single layer) was computed using the plane-wave DFT code ABINIT [33,154]. The same pseudopotentials, exchange correlation functional and electronic temperature were used between ABINIT and Helical DFT. Energy cutoffs between 40 and 60 Ha, 30 × 30 × 1 k-points, and a cell vacuum of 25 Bohr in the direction orthogonal to the sheets, were employed. At the end of the geometry relaxation procedure, the atomic forces and the cell stress were of the order of 10 −5 Ha/Bohr and 10 −8 Ha/Bohr 3 , respectively. The agreement of these parameters with existing literature is quite good [47], thus lending confidence to the physical properties of the X nanotubes as revealed via our simulations. We apply a twist to each of these systems by setting α between 0.003 and 0.006 (this corresponded to between 2.47 and 8.86 degrees/nm of imposed rate of twist). With all the other parameters of the computational method fixed to values described earlier, the only remaining quantities that can dictate the accuracy of the numerical solutions are fineness of the real and reciprocal space discretizations. Accordingly, we study the convergence behavior of the ground state energy and the atomic forces as a function of the mesh spacing h, and the number of reciprocal space points N η used in the calculations. The results are shown in Figure 5. For the mesh convergence study, we used h = 0.15 Bohr to evaluate the reference value while computing errors, while for studies involving convergence with respect to reciprocal space discretization, we used N η = 21 as reference.
From the figures, we see that the numerical method converges systematically in each of the cases under study. By fitting straight lines to the convergence data with respect to h, we observed slopes between 5.5 and 6.5 which are somewhat lower than values observed for finite difference calculations using (untwisted) cylindrical coordinates [47]. We are also able to estimate that a mesh spacing of about h = 0.3 Bohr, and a value of N η = 15 are more than sufficient to reach chemical accuracy thresholds in all cases (i.e., 10 −3 Ha/atom in the energies and 10 −3 Ha/Bohr in the atomic forces), and we used these discretization choices in structural relaxation calculations in subsequent sections. Figure  6 shows the consistency of the forces and the energies as computed by Helical DFT at this level of discretization (i.e., numerical derivatives of the free energy per unit cell as computed via eq. 29, yield the atomic force as computed via eq. 61). To compute the energies and band structures of relaxed structures, we employed the finest discretization parameters that we could reliably afford within computational resource constraints. This corresponded to the choices h = 0.25 Bohr and N η = 21.
Next, we come to a discussion of verification of the numerical method against results produced by standard, widely used plane-wave codes such as ABINIT [33,154]. As described earlier, this can be an arduous endeavor since such codes may require a very large number of atoms to be included in the periodic unit cell, in order to mimic the systems being simulated via Helical DFT. Moreover, in order to accurately accommodate the boundary conditions implemented in Helical DFT, a large amount of vacuum padding has to be often employed in the plane-wave code unit cell, and nanotube-like structures tend to encase a large amount of vacuum as it is. These factors together can result in slow convergence of the electrostatics problem, as well as, poor conditioning of the systems of equations being solved by the plane-wave code. The latter issue, in turn, leads to SCF convergence problems which tend to worsen if calculations at high accuracies are required (i.e., upon using a large value of E cut for the plane-wave code). With these factors in mind, we chose the armchair carbon and silicon nanotube systems described above for comparison against ABINIT. For the former, we did not prescribe any twist and use a 64  atom unit cell. For the latter, we prescribed a twist of α = 0.1, and used a 360 atom unit cell. While dealing with these systems in ABINIT, periodicity was naturally enforced along the Z axis, Dirichlet boundary conditions were enforced along the X and Y axes by padding with a large amount of vacuum, and an SCF preconditioner (diemac option in ABINIT) was used to deal with instabilities associated with spatial inhomogeneities in the periodic unit cell. Helical DFT was made to use four atom unit cells for both examples.
For each of these model systems, we observed that the energies (in Ha/atom) and the forces (in Ha/Bohr), from ABINIT and Helical DFT agreed with each other to O(10 −4 ), thus giving us confidence in the accuracy of the results produced by our method. 10 Based on the above tests, we were also able to observe that even a well optimized plane-wave code like ABINIT can take up to orders of magnitude more in simulation time (measured in c.p.u. hours) compared to Helical DFT, when simulations of nanotube structures (particularly, ones with imposed twist) are desired. This makes our computational method a particularly attractive choice in the first principles characterization of such systems. The relative efficiency of our method stems from the use of a coordinate system and a computational domain that are well adapted to the geometry of the twisted structure, and also from the appropriate use of symmetry. To highlight the latter aspect, we considered again the silicon nanotube system subjected to a twist of α = 0.1. We used Helical DFT to calculate the ground state electronic structure of this system by considering the following four equivalent scenarios: The single core wall times required for each SCF step, and computation of the atomic forces at the end of the SCF iterations are compared in Figure 7. From these plots, it is clear that the SCF wall time is approximately 20 times lower for the case with full symmetry adaptation, when compared to the case in which no cyclic or helical symmetries were used. Even more drastically, the computational wall time for the calculation of the force is about 3 orders of magnitude lower for the former case, when compared to the latter. These computational advantages tend to be even more dramatic for simulations in which the angle of twist is relatively low (e.g. α = 0.0005 to 0.005), and such cases tend to arise routinely while probing the torsional response of the nanotubes in the linear elastic regime, as described in the next section.
Finally, we show in Figure 8 the strong scaling behavior of the numerical implementation. We use case (d) described above for this study. We see that up to 16 computational cores, the code has a strong scaling efficiency of about 60 %. This follows the strong scaling efficiency of the CheFSI step closely, since this forms the dominant computational cost in every SCF step (see Figure 7(a)). The scaling of the force computation step is far worse, dropping to about 10 % at 16 cores. In general, the scaling behavior is expected to improve for problems with a larger number of η and ν points (e.g. for simulations of nanotubes of large diameter) since the current version of the code only uses parallelization over different values of η and ν. Improvement of the scaling behavior of the code, particularly by use of domain decomposition and band parallelization techniques in conjunction with the MATLAB Parallel Server framework (to enable deployment over distributed memory computers) is the scope of future work. 11

Computation of torsional stiffness from first principles
We now turn to demonstrations of the use of our computational method for evaluation of materials properties from first principles. We first concentrate on the mechanical response and evaluate the torsional stiffness of the X nanotubes in the linear elastic regime, ab initio. We choose 9 to 10 nanotubes of each material, about half of which are of zigzag type and the other half armchair. The nanotubes all had radii in the range 1 to 3 nm, approximately. To carry out these simulations, we choose a four atom unit  cell for the untwisted nanotube in each case, and perform structural relaxation using the FIRE algorithm [98] till all force components on all the atoms in the simulation cell dropped below 10 −3 Ha/Bohr. We then successively increase α to impose twist, and in each case re-perform structural relaxation (see Figure 9 for some examples of the relaxation procedure). To avoid the appearance of torsional instabilities, we ensured that  the prescribed rate of twist on the system was less than about 4.5 degrees per nanometer [36], and this corresponded to choosing α between 0.0005 and 0.005. We express the amount of applied twist per unit length of the tube (i.e., the rate of twist) as β = 2πα τ , and compute the twisting energy per unit length of the structure as the difference in the ground state free energy per unit fundamental domain between the twisted and untwisted configurations (after atomic relaxation is carried out in both cases), i.e.: Here, G| β and G| β=0 denote the symmetry groups associated with the twisted and untwisted structures, respectively. Also, P * * and P * denote the collections of positions of the atoms in the fundamental domain, after relaxation in each case. For each of the nanotubes, we verified that mechanical response was in the linear regime, by fitting U twist (β) to a function of the form U twist (β) = c × β q and observing that q ≈ 2.0 holds. We show a few examples in Figure 10. Next, using the above data, we estimated the twisting stiffness of each nanotube, defined as: For each category of nanotube (i.e., armchair or zigzag, and type of material), we then studied the variation of k twist with the nanotube radius (computed as the average of the radial coordinates of all atoms in the fundamental domain), by using a fit of the form: The results from this procedure are shown in Figure 11 and the values of κ and p obtained in each case are displayed in Table 2. Note that generation of this torsional response data required hundreds of individual simulations, which would not have been possible without the use of a specialized computational method such as the one presented here. A few comments are in order at this stage. First, we observe that the value of the exponent p is nearly 3 in every case. This suggests that the torsional response of the   Table 2: Torsional stiffness parameters for the X nanotubes, with k twist = κ × R p tube tubes is consistent with linear elasticity theory, in which k twist for a thin elastic tube with unit length, radius R tube , thickness t, and shear modulus G can be expressed as GtπR 3 tube [158]. From this, it is possible to estimate the thickness-normalized shear modulus (i.e., Gt) of the Xene sheets as κ/π. Second, by comparing the different values of κ, we see that they span an order of magnitude across the different elements. In particular, for a given radius, k twist is the highest for carbon nanotubes and the lowest for those of tin, while nanotubes of silicon and germanium have intermediate values of this quantity close to each other. Third, for each material, the torsional response is quite similar in the armchair and zigzag directions with variations less than about 1.5%, except for the case of tin, in which case the variation is more substantial. This largely isotropic torsional response for the Xene nanotubes is quite distinct from the bending response of their sheet counterparts, which show strong anisotropic behavior that is correlated with the value of the normalized buckling parameter (i.e., δ/a) for each material [47]. Our findings on the mechanical response of carbon nanotubes under torsion are broadly consistent with earlier studies for this material that used empirical potentials or tight-binding calculations [36,38], although the value of κ reported here is lower from [36], where Tersoff potentials were used [159]. Finally, we mention in passing, the effects of atomic relaxation. In general, if relaxation is not performed after the imposition of twist, the value of k twist for the system tends to be higher. The degree of variation can be quite different depending on the material involved. For carbon nanotube systems, we observed that k twist for an unrelaxed system was typically higher by a factor of about 1.08, whereas for silicon nanotubes, this factor had the higher value of about 1.38. Generally, these higher values of k twist also imply higher values of κ by the same factors, although the value of the exponent p continues to be about 3, when the fitting in eq. 84 is used.

Investigation of electronic properties of nanotubes undergoing torsional deformation
We now discuss the variation in electronic properties of nanotubes as they are subject to twisting. Due to the ability of our computational method to use symmetries connected with the system, electronic band-diagrams along both η and ν can be obtained from Helical DFT. Moreover, the eigenvalues λ j (η, ν) as j is held constant and η, ν are varied, can be plotted as a two-dimensional surface. Since η and ν serve to label the set of characters, and are natural quantum numbers for twisted structures, they serve to provide a clean and intuitive interpretation of the electronic states of the system, and allow easy identification of the size and type of band-gaps. In contrast, the traditional band diagram for a quasi-one-dimensional system using a periodic method can be far more complicated, even for an untwisted structure. We show some examples of this contrast in Figures 12  and 13. Armed with the above tools, we study the variation in the bandgaps of nanotubes as they are subjected to twisting. For reasons explained later, we mainly concentrate on investigations related to armchair X nanotubes, although we also briefly comment on our findings related to zigzag X nanotubes subsequently. The behavior of carbon armchair nanotubes in particular, has received much attention in the literature [62,135,143,144], and serves as an important benchmark against which our results can be validated. Such nanotubes are known to be metallic [143,144,160] although in practical calculations, a vanishingly small bandgap at the location η = 1 3 , ν = 0 (or equivalently, η = − 1 3 , ν = 0)

34
(a) 2D surface plot of the eigenvalues λ j (η, ν), for j = 8.  may be observed [47]. Upon twisting, armchair carbon nanotubes undergo a metal-tosemiconductor transition, with the characteristic feature that the bandgap-versus-rateof-twist plot has a slope of 3 t 0 R tube in the linear regime (i.e., in the neighborhood of zero twist). Here t 0 is the tight-binding hopping parameter for carbon [62]. Using armchair carbon nanotubes of radii 1.08, 1.48 and 1.88 nm as examples, we used Helical DFT to compute the slope of the bandgap-versus-rate-of-twist plot in the linear regime and obtained values of t 0 between 2.6 and 3.0 eV (see Figure 14). These agree well with the literature [62,144,161], giving us confidence in the quality of our subsequent simulations. Upon twisting these nanotubes further, the band gap is known to further increase and then decrease, as the tube alternates between metallic and semiconducting states, and the period of oscillation (of the band gap versus rate of twist plot) is theoretically known to be [62,135,143,144]: Here a denotes the carbon-carbon bond length (see Table 1). Using Helical DFT, we were able to compute the electronic density of states near the Fermi level and qualitatively verify the metal-to-semiconductor transitions in the armchair carbon nanotubes as they are twisted (see Figure 16(a)). To verify that Helical DFT also reproduces the quantitative aspects of the variation, we fit the band gap data from Helical DFT, to a general sine curve of the form: band gap = s 1 sin 2πα from which, the period of oscillation may be computed as: We verified that ξ fit period and ξ theory period are in close agreement in all cases under study (see Figure 14 for a specific example). An alternate means of quantifying this agreement, following [62], is to equate ξ fit period and ξ theory period , and estimate the bond length a, from this instead. In other words, by writing: or more generally, we may evaluate the exponent µ and the constant σ from a plot of s 2 versus R tube , and from this, we may further estimate the bond length as: Using this procedure, we arrived at µ = −1.98, and a fit = 1.37 angstrom, both of which are very close to the expected values of −2.00 and 1.40 angstrom, respectively. These results give us further confidence in the quantitative results obtained using Helical DFT.
Turning to the broader class of armchair group IV nanotubes (i.e., X = Si, Ge, Sn) we make the following observations using the data obtained from Helical DFT. In general, these nanotubes are semiconducting, with a direct band gap located at the same position as the armchair carbon nanotubes, i.e., η = 1 3 , ν = 0 (or equivalently, η = − 1 3 , ν = 0) for untwisted tubes. Upon twisting, these tubes also undergo periodic oscillations in their band gaps, 12 although the amplitudes of the oscillations are generally more muted than 12 The location of the band gap initially continues to be the same as that of the untwisted tube, but then it transitions to small values in ν (i.e. ν = 1, 2, etc.), while remaining at the same location in η (i.e., η = 1 3 ). Thus, for relatively small twists, the band gap continues to be a direct one. Upon further application of twist however, the band gap becomes indirect and the eigenvalue just above the Fermi level is associated with a different value of ν as compared to the eigenvalue just below the Fermi level, although the value of η associated with these eigenvalues continues to be 1 3 . Data points computed using Helical DFT Sine-curve fit Linear fit near zero twist Figure 14: Analysis of the variation of band-gap with applied twist, using an armchair carbon nanotube example (Radius = 1.07 nm). The straight line fit near zero enables the evaluation of the tight-binding hopping parameter t 0 , which comes out to be 2.897 eV, in close agreement with [62,144,161]. The sine curve fit (in the non-linear response region) enables evaluation of the periodicity in the band gap variation and yields ξ fit period = 0.1154 rad/nm. The theoretical value from eq. 85 is ξ theory period = 0.1217 rad/nm, in close agreement.
the case of armchair carbon nanotubes, and we did not observe metal-to-semiconductor type transitions for most tubes. For tubes with larger radii however, the untwisted states can be associated with vanishingly small band gaps to begin with -owing to the decay relations obeyed by the band gaps [47,132], and these tubes are likely to be practically metallic at room temperature. Therefore, changes to the band gap upon application of twist can be more easily discerned (See Figures 13 and 15 for an example involving an armchair silicon nanotube). To quantify the periodic changes in the band gaps, we obtained the period of oscillation in each case using the sine curve fitting procedure outlined above (eq. 86), and computed the power law dependence of the period on the tube radius by means of eq. 89 (see Figure 17). The values of c and µ so obtained are shown in Table 3.
The results are clearly suggestive of the fact that the period of variation of the band gap scales in an inverse quadratic manner with the tube radius for all armchair X nanotubes. We also observed that evaluation of eq. 90 using the values of σ shown in Table  3 leads to quantities that are fairly close to the values of a shown in Table 1, for each armchair X nanotube, suggesting that the theoretical relation in eq. 85 is generally valid for this entire class of nanotubes.
Finally, we touch upon our investigations related to zigzag X nanotubes. These can be of different "types" [47,143], i.e., Type I, II or III, depending on whether mod(N, 3) = 1, 2 or 0. In general, zigzag X nanotubes, barring Type III carbon variants, are semiconducting [47,162], and the untwisted tubes have direct bandgaps located at the following values of η and ν -Type I carbon nanotubes: η = 0, ν = N−1 3 ; other Type I nanotubes: η = 0, ν = N+2 3 ; Type II nanotubes: η = 0, ν = N+1 3 ; Type III nanotubes: η = 0, ν = N 3 . We found that the band gaps of Type I and II zigzag X nanotubes tend to have a rather limited response to torsional deformations, consistent with earlier observations made re-  Table 3: Parameters for the scaling law s 2 = σ ×R µ tube for armchair X nanotubes. Here, s 2 is the bandgap oscillation parameter as defined in eq. 86. The value of µ in each case is close to −2.00, suggesting that the period of variation of the band gap scales in an inverse quadratic manner with the nanotube radius for these tubes.
(a) 2D surface plot of the eigenvalues λ j (η, ν), for j = 8.   garding zigzag carbon nanotubes specifically [135,143,144]. For most of these types of materials, the band gaps are non vanishing at zero twist for even relatively large radii tubes and the subsequent changes to their band gaps due to twisting are fairly small at the levels of torsional deformation we considered. This tends to cause issues in discriminating between actual changes to the band gaps due to deformation, and the numerical noise associated with the simulations. Therefore, although we did observe oscillatory patterns in the band gap versus rate of twist plots (see Figure 18 for an example) we found it difficult to extract scaling laws from this data unambiguously. Out of all the different zigzag X nanotubes however, the Type III variants of carbon are metallic, especially at larger radii (i.e., when curvature effects are minimal) [47,163], and we observed such tubes to be quite sensitive to torsional deformations. Similar to the case of armchair nanotubes, we observed these tubes to show oscillatory behavior between metallic and semiconducting states (see Figure 18), and an analysis of the period of variation of the band gap (using eq. 86 and 89) yielded µ = −1.98, thus suggesting an inverse quadratic dependence on the radius. A thorough re-investigation of scaling laws in the electronic response of zigzag X nanotubes, by making use of more accurate numerical techniques (based on spectral methods [164,165], for instance) remains the scope of future work.

Conclusions
In summary, we have presented a computational technique that allows systems associated twisted geometries to be simulated efficiently and accurately from first principles. We have formulated the symmetry adapted governing equations, laid out numerical implementation strategies and detailed various aspects of our implementation. Our technique uses a higher order finite difference discretization scheme based on helical coordinates, employs ab initio pseudoptentials and can be used to simulate quasi-one-dimensional systems, as well as their deformations, conveniently and without needing major computational resources. As an application of our method, we have systematically studied  ) and (c) include data from Helical DFT, as well as sine curve fits (dotted lines) used to determine the band gap oscillation parameter s 2 (eq. 86). Sub-figure (d) explores the variation of this parameter with the tube radius (eq. 89). The slope of each of the straight line fits is close to −2.00, suggesting that the period of variation of the band gap scales in an inverse quadratic manner with the nanotube radius.
the behavior of single wall zigzag and armchair group-IV nanotubes in the range of (approximately) 1 to 3 nm radius, as they undergo twisting. Through an extensive series of simulations, we have demonstrated how certain mechanical properties of these nanotubes can be extracted from first principles using our technique, and we have also elucidated different aspects of the variation in the electronic properties of these materials as they undergo torsional deformation. In particular, using our simulations, we have been able to extend some well-known features of the electro-mechanical properties of carbon nanotubes to the broader class of Group IV nanotubes.
As a follow up of this work, we aim to employ the computational technique discussed here for the study of other nanotube materials, including multi-wall elemental nanotubes, and those made from transition metal dichalcogenides. An efficient C/C++ implementation of the computational method which makes use of domain decomposition and band parallelization (in addition to the currently implemented parallelization in η and ν), to improve scaling and computational wall time performance is the scope of ongoing and future work. Concurrently, the development of an efficient spectral scheme [164,165] in the spirit of [166], which overcomes some of the inherent limitations of the current finite difference technique is also an area of active investigation. Finally, a long term goal associated with applications of the current computational method involves the design and discovery of exotic materials phases which show strong coupling between mechanical deformations (such as twist and extension/compression) and other electronic/optical/magnetic/transport properties.