Ab initio framework for systems with helical symmetry: theory, numerical implementation and applications to torsional deformations in nanostructures
AAb initio framework for systems with helical symmetry: theory,numerical implementation and applications to torsionaldeformations in nanostructures
Amartya S. Banerjee a, ∗ a Department of Materials Science and Engineering, University of California, Los Angeles, CA 90095,U.S.A
Abstract
We formulate and implement Helical DFT — a self-consistent first principles simulationmethod for nanostructures with helical symmetries. The mathematical framework forclassifying nanostructures shows that for a vast class of nanomaterials — including im-portant technological materials such as nanotubes, nanosprings and nanowires, as well asmiscellaneous examples from chemistry and biology — the spatial atomic arrangementpossesses helical symmetries. The overwhelming preponderance of such helical structuresin all of science and engineering and the likelihood of these systems being associatedwith exotic materials properties makes them one of the most important and interestingcategories of materials. This provides the motivation for the development of accurate andpredictive simulation techniques for their study.The electronic states in a nanostructure with helical symmetries can be characterizedby means of special solutions to the single electron problem called helical Bloch waves .We use arguments from Fourier analysis and the theory of elliptic operators to rigorouslydemonstrate the existence and completeness of such solutions. We describe how theKohn-Sham Density Functional Theory (KS-DFT) problem for a helical nanostructurecan be reduced to a fundamental domain with the aid of helical Bloch waves. We thensystematically derive the governing equations. A key component in our mathematicaltreatment is the definition and use of a helical Bloch-Floquet transform to perform a“block-diagonalization” of the Hamiltonian in the sense of direct integrals. We develop asymmetry-adapted finite-difference strategy in helical coordinates to discretize the gov-erning equations as posed on the fundamental domain, and obtain a working realizationof the proposed approach. We verify the accuracy and convergence properties of ournumerical implementation through examples.Finally, we employ Helical DFT to study the properties of zigzag and chiral single wallblack phosphorus (i.e., phosphorene) nanotubes. We use our simulations to evaluate thetorsional stiffness of a zigzag nanotube ab initio . Additionally, we observe an insulator-to-metal-like transition in the electronic properties of this nanotube as it is subjectedto twisting. We also find that a similar transition can be effected in chiral phosphorenenanotubes by means of axial strains. The strong dependence of the band gap of thesematerials on various modes of strain suggests their possible use as nanomaterials withtunable electronic and transport properties. Notably, self-consistent ab initio simulationsof this nature are unprecedented and well outside the scope of any other systematic firstprinciples method in existence. We end with a discussion on various future avenues andapplications. 1 a r X i v : . [ phy s i c s . c o m p - ph ] A ug eywords: Kohn-Sham density functional theory, helical symmetry, phosphorene,nanotube, torsional deformations.
1. Introduction
The discovery and characterization of novel nanomaterials and nanostructures con-stitutes one of the principal areas of scientific research today [1, 2]. Such materials andstructures hold the promise of unlocking remarkable and unprecedented material prop-erties that are otherwise unavailable in the bulk phase (i.e., crystalline materials). Inrecent years, the discovery of novel nanostructures has garnered much attention and ac-claim [3, 4], and the unusual properties of these new materials have led to ground breakingapplications in almost every branch of science and engineering [5, 6].Nanostructures appear in various morphologies (including fullerenes, nanotubes, nan-oclusters and two-dimensional materials), and are usually associated with non-periodicsymmetries. The mathematical framework for classifying nanostructures [7, 9, 10] showsthat a vast class of these materials can be described as being helical , i.e., their spatialatomic arrangement possesses helical symmetries. Helical structures include importanttechnological materials such as nanotubes (of any chirality), nanoribbons, nanowires andnanosprings; miscellaneous chiral structures encountered in chemistry; and examples frombiology, including tail sheaths of viruses and many common proteins [7, 11]. Figure 1shows instances of helical structures that have been actively investigated in the literature.Helical structures have been conjectured to be a fertile source of novel materialswith unusual and attractive properties [7]. This is due to the fact that atoms in suchstructures find themselves in locally similar environments [7]. Coupled with the quasi-one-dimensional nature of these systems, as well as the presence of symmetries in theunderlying governing equations, this makes it likely that collective or correlated electroniceffects (such as those leading to ferromagnetism, ferroelectricity and superconductivity)can emerge in these materials [15]. Contrarily, helical structures are also inherently chiraland can therefore serve as natural examples of materials systems in which certain formsof symmetry breaking in the governing equations can lead to unconventional transportphenomena [16–18].Given the relative abundance of helical nanostructures in existing materials, their like-lihood of being associated with hitherto undiscovered forms of matter displaying exoticmaterials properties, and their overall scientific and technological importance, there ap-pears to be a pressing need for reliable and efficient computational tools for studying suchsystems. The broad goal of the present contribution is to take important foundationalsteps in addressing the above scientific issue. Specifically, we present here the mathemat-ical formulation and numerical implementation of a novel computational method called ∗ Corresponding author
Email address: [email protected] (Amartya S. Banerjee) Atomistic and molecular structures with non-periodic symmetries have been termed as
ObjectiveStructures in the mechanics literature [7]. First principles calculations for such structures was the topicof investigation of [8] and the current contribution continues and extends that line of work, i.e., it canbe viewed as a particular flavor of
Objective Density Functional Theory . Preprint submitted to Elsevier August 6, 2020 igzag Carbon Nanotube Armchair Carbon Nanotube (a) Nanotubes of different chirality.
Silicene Nanoribbon (Flat) Phosphorene Nanoribbon (Twisted) (b) Nanoribbons od different 2D materials.
Amyloid Protein Fibrils (from di ff erent animal and plant sources) (c) Molecules of biological origin.Figure 1: Examples of helical structures from nanotechnology and biology. Images of the Amyloid proteinwere obtained courtesy of The Protein Data Bank [12–14]. Helical DFT, that can simulate helical structures ab initio . We also obtain a practicalworking realization of this (density functional theory based) self-consistent first princi-ples technique, and illustrate some of its capabilities through the study of an emergentnanotube material with interesting properties.To put our work into perspective, we remark that the use of first principles (i.e.,quantum mechanical) techniques to design and study materials is a very active areaof scientific endeavor today, and it forms the bulk of computational materials scienceresearch [19–23]. Among the wide array of first principles methods available, Kohn-Sham Density Functional Theory (KS-DFT) [24] enjoys widespread usage since it offersa good balance between computational cost and physical accuracy as compared to othertechniques [25]. The pseudopotential plane-wave method, also called
Plane-wave DFT ,is the most widely used implementation of Kohn-Sham theory [26–29], and it involvesexpanding the unknowns into linear combinations of plane-waves. Since plane-wavesare naturally associated with periodic symmetries (they are in fact eigenfunctions oftranslational symmetry operators), Plane-wave DFT is ideally suited for studying bulk(i.e. periodic or crystalline) systems, and is often found to be fundamentally inadequate3or studying systems with non-periodic symmetries. In particular, using a Plane-waveDFT code for studying a helical structure such as a chiral nanotube can require the useof large periodic unit cells often containing many hundreds (or even thousands) of atoms. In contrast, a computational method which is faithful to the underlying helical symmetryof such a structure would require a small helical unit cell, containing far fewer numberof atoms. Since ground state electronic structure calculations using density functionaltheory (DFT) scale as the cube of the number of atoms in the unit cell, while excitedstate calculations scale as the fourth power, the difference in simulation run times forsuch calculations, in these two scenarios (i.e., the correct use of helical symmetry vs.incorrect use of periodic symmetry) can be drastically different in practice.The above considerations form our point of departure from a conventional formulationand implementation of KS-DFT, to one that is adapted for helical systems. In order toformulate the equations of KS-DFT for a helical unit cell, an appropriate version of theBloch Theorem [37, 38] is required. We establish this result rigorously in this work, anduse it to set up an electronic band theory for helical structures. Subsequently, we developthe notion of helical Bloch states, and use their properties to derive of the equations ofKS-DFT, as they apply to helical systems. Our use of rigorous mathematical argumentsand appropriate mathematical tools is one of the highlights of our framework, and itallows the governing equations to be obtained systematically, and without recourse toan excessive amount of intuition. As far as we are aware, our work is the first inpresenting such a derivation, and also in expressing the detailed form of the equations ofKohn-Sham theory for helical structures. The final form of the equations are such thatthey are readily suited for implementation within systematically convergent electronicstructure methods such as those based on finite differences [30–33], finite elements [34–36] or spectral basis functions [8, 44, 45]. We choose a symmetry adapted finite differencemethod in helical coordinates for discretizing the governing equations in this work, andset up a computational framework for numerically solving the discretized equations in aself-consistent manner. This gives us a working realization of an ab initio computationaltool — called Helical DFT — that can be used to perform predictive simulations ofhelical systems in a systematic and efficient manner. It can therefore aid in the discovery,synthesis and characterization of helical structures. Subsequently, the remainder of thiswork focuses on illustrating various numerical and application oriented aspects of thisnovel computational tool through examples based on nanotube systems. To the best ofour knowledge, Helical DFT is the first computational method for helical systems that In contrast to plane-waves, the use of real space techniques based on finite differences [30–33] orfinite elements [34–36] allow for non-periodic boundary conditions to be imposed in a straight-forwardmanner. However there does not appear to be any prior work on using these techniques for self-consistentfirst principles calculations of helical systems. Due to the infinite nature of helical groups, the mathematical arguments presented here are ofsomewhat different and more subtle nature as compared to the ones that can be employed for cyclicgroups [39]. However, they can be seen as being broadly connected in the sense that they both deal withFourier analysis of the respective symmetry groups [40, 41]. Since a rigorous thermodynamic limit theory for the Kohn-Sham problem is unknown [42, 43], aderivation of the equations of the theory, as it applies to condensed matter systems often makes use ofphysical intuition. This process is prone to conceptual errors however, and we are aware of literaturethat lists certain terms of the equations incorrectly. In any case, the final form of the equations appearto be well known in the electronic structure community at large, since DFT codes routinely make use ofthem for simulating the crystalline phase.
4s based on first principles, and one that also behaves systematically with respect toconvergence properties. This, among other reasons, is made possible by our use of theaforementioned helical coordinate system. To the best of our knowledge, this has notbeen employed in electronic structure calculations before.While the study of helical structures has much scientific and technological merit in ofitself, the development of a computational method for studying such systems also bringswith it the added benefit of being able to simulate the behavior of nanomaterials undertorsional deformations. As explained in [7, 39], homogeneous deformation modes arecommensurate with periodic symmetries (i.e., applying a homogeneous deformation toa periodic structure results in another periodic structure), while certain inhomogeneousdeformation modes can be associated with non-periodic symmetries. An attempt tostudy such inhomogeneous deformations while using a periodic method is likely to involvevarious uncontrolled approximations, complications and computational inefficiencies [46,47]. This issue appears to have been recognized for some time in the nanomechanics andmaterials literature, leading to a considerable body of work centered around suggestionspresented in [7], whereby pure bending deformations in atomistic systems are simulatedusing cyclic symmetries, while helical symmetries are used to simulate torsion [48–55]. Apersistent issue with the simulations in these studies however, is that they have all beencarried out using interatomic potentials or tight binding methods. Due to the well knowndeficiencies of these techniques in simulating real materials [47, 56–59], true first principlessimulation methods that behave systematically, and also take into account cyclic and/orhelical symmetries have been deemed highly desirable [7, 8, 55]. There has been recentprogress on this very issue with regard to cyclic symmetries [39, 60], and the resultingcomputational methods have been used to study the bending behavior of nanoribbonsand sheets of two dimensional materials ab initio . In this sense, the current contributionfollows up on this line of work by making a first principles simulation framework fortorsional deformations available. Consequently, through the use of this framework, weare able to extract the behavior of nanotubes of black phosphorus (i.e., phosphorenenanotubes) and study their mechanical and electronic response as they are subjected totwisting. The coupling of these responses leads to some interesting electronic transitionsin this material that is likely to make it an attractive candidate for sensing, modulationand actuation applications.The rest of this work is organized as follows. Section 2 establishes the mathematicalframework for a systematic formulation of the governing equations, and also derivesthe relevant expressions explicitly. Section 3 discusses formulation of a numerical schemebased on finite differences in helical coordinates, and Section 4 presents simulation studies.Section 5 summarizes the work and suggests avenues for future research. The appendices Exploitation of helical symmetries in ab initio calculations has also been considered in the chemistryliterature in the context of Linear Combination of Atomic Orbitals (LCAO) methods [61–66]. However,these methods differ in their perspective from the current contribution in that they concentrate on usingsymmetry-adapted basis functions for reducing the computational cost of the multi-center integrals andthe Hamiltonian matrix elements, whereas our focus is on the formulation of symmetry-adapted cellproblems (in helical coordinates), and a systematically convergent numerical treatment of these cellproblems. Due to basis incompleteness and superposition errors, it is often non-trivial to systematicallyimprove the quality of the numerical solutions obtained via LCAO methods, in contrast to the techniquespresented here. Finally, the connection of helical symmetries with torsional deformations, as well as theeffect of such deformations on other material properties does not appear to have been considered in thechemistry literature.
2. Formulation
In this section, we describe the key aspects of Helical DFT. We begin with a formaldiscussion of helical groups and helical structures in Section 2.1, and then discuss Kohn-Sham DFT, as it applies to such systems in Section 2.2. The atomic unit system with m e = 1 , e = 1 , (cid:126) = 1 , π(cid:15) = 1 , is chosen for the rest of the work, unless otherwisementioned. A helical structure (i.e. a structure with helical symmetries) can be defined throughthe action of a helical group on a set of non-degenerate points in space. This definitionmakes it necessary for us to make the notion of a helical group precise. Following standardpractice in the literature [7–10, 48, 67, 68], we introduce helical groups as subgroups ofthe Euclidean group in three dimensions. This requires us to introduce some relevantnotation and basic rules regarding operations with isometries, as we now do.
Let e , e , e denote the standard orthonormal basis of R and let ( x , x , x ) denotethe Cartesian coordinates of a generic point x ∈ R . An isometry (or rigid body motion)in R will be denoted using the notation Υ = ( R | c ) , with R ∈ SO (3) denoting therotation part of the rigid body motion, and c ∈ R denoting the translation part. Theaction of Υ : R → R on a point x ∈ R is written as Υ ◦ x = Rx + c . Given a collectionof points S ⊂ R , we will use the notation Υ ◦ S to denote the action of the isometry oneach of the points in S , i.e., Υ ◦ S := (cid:91) x ∈ S Υ ◦ x . (1)There is a natural multiplicative operation associated with isometries (denoted as • here) that arises as a composition of their maps. Specifically, given isometries Υ =( R | c ) , Υ = ( R | c ) , we may define a third isometry Υ = Υ • Υ such that Υ ◦ x =(Υ • Υ ) ◦ x := Υ ◦ (Υ ◦ x ) . It follows that Υ = ( R R | R c + c ) , and that in generalthe operation • is not commutative (due to non-commutativity of finite rotations aboutarbitrary axes). The • operation also allows the definition of whole number powers of Υ ,i.e., for n = 1 , , . . . , we may define Υ n := Υ • Υ • Υ . . . ( n times ) . It is then easy to check We will use the following notation in what follows: f ( · ) will be used to denote a function when wedo not wish to highlight the dependence of the function on its arguments. (cid:107)·(cid:107) will be used to denote thenorm of a function or vector and (cid:104)· , ·(cid:105) will be used to denote the inner product. Often, we will attacha subscript to these symbols to denote the specific space in which the norm or inner product is beingconsidered. Vectors and matrices in R or R will be denoted in boldface, with lower case letters reservedfor vectors and uppercase letters used for matrices. We will sometimes use the · symbol between vectorsin R or R , to denote the inner product. As the name suggests, isometries preserve distances (and hence, also angles), i.e., ∀ x , y ∈ R , and ageneric isometry Υ , it holds that (cid:107) Υ( x ) − Υ( y ) (cid:107) R = (cid:107) x − y (cid:107) R . Υ n admits the expression Υ n = (cid:0) R n (cid:12)(cid:12) ( n − (cid:88) j =0 R j ) c (cid:1) , where the notation R is used todenote the identity matrix.The identity isometry leaves every x ∈ R invariant and can be written as ( I | ) ,with I denoting the identity matrix and denoting the null vector in R . Given theisometry Υ = ( R | c ) , we can form the isometry Υ (cid:48) = ( R − | − R − c ) , which satisfies Υ • Υ (cid:48) = Υ (cid:48) • Υ = ( I | ) . Hence, we will denote Υ (cid:48) as Υ − — i.e., the inverse isometryto Υ . The set of all isometries so defined, i.e., E = { Υ = ( R | c ) : R ∈ SO (3) , c ∈ R } ,together with the operation • and the inverse element defined above, form a group [69]. Let α and τ be real numbers such that ≤ α < and τ > , and let R πα denote arotation around axis e by angle πα . Then, the rigid body motion Υ h = ( R πα | τ e ) willbe called a helical isometry about axis e . The action of Υ h on a point x ∈ R is torotate it by angle πα about axis e , while also translating it by τ along the same axis. Furthermore, applying the formulae for the powers of isometries and their inverses shownabove, we see that for m = 1 , , . . . , Υ m h = ( R πmα | mτ e ) and Υ − h = ( R − πα | − τ e ) .Combining these, we may define Υ m h for any m ∈ Z as Υ m h = ( R πmα | mτ e ) , with the m = 0 case automatically resulting in the identity isometry ( I | ) . We may thereforestate: Proposition 2.1 ( Helical group generated by a single element ) . The set of isome-tries G = (cid:8) Υ m h = ( R πmα | mτ e ) : m ∈ Z (cid:9) , (2) forms a group under the operation • . Additionally, let N ∈ N , let Θ = 2 π N and for n = 0 , , . . . , N − , let R n Θ denotea rotation around axis e by angle n Θ . Then the set of isometries endowed with theoperation • C = (cid:8) Υ n c = ( R n Θ | ) : n = 0 , , . . . , N − (cid:9) , (3)forms a cyclic group [39] of order N . Note that since the rotational parts of the isometriesin group G and C all share e as the common axis of rotation, the elements of G and C commute (i.e., for any Υ m h ∈ G and Υ n c ∈ C , Υ n c • Υ m h = Υ m h • Υ n c holds.) . We may nowconsider the direct product of the groups G and C defined above to obtain a new helicalgroup : Since only pure rotations are included, this is the so called Euclidean group of direct isometries inthree dimensions [69]. The full Euclidean group also includes improper rotations. Most of the discussion in this work naturally also extends to the case when − < α < . However,we will not be considering that case here. Alternately referred to as a screw transformation in the crystallography literature [9]. A simple way to see this is to resolve x along and perpendicular to e , i.e., x = x e + x ⊥ e ⊥ , where (cid:104) e , e ⊥ (cid:105) R = 0 and (cid:107) e ⊥ (cid:107) R = 1 . Then, Υ h ◦ x = ( x + τ ) e + x ⊥ ( R πα e ⊥ ) . While the discussion presented here already makes it evident that both G and G are groups (andare in fact Abelian groups), see [10] for a more complete derivation of these groups, as well as other typesof helical groups not considered in this work. Note that the groups G and G contain a group of translations as a normal subgroup if α is a roposition 2.2 ( Helical group generated by two elements ) . The set of isometries G = (cid:8) Υ m h • Υ n c = ( R πmα + n Θ | mτ e ) : m ∈ Z ; n = 0 , , . . . , N − (cid:9) , (4) forms a group under the operation • . Since G and C are generated by single elements, they are Abelian groups. Further-more, since G is generated by two elements (i.e., the generators of G and C ) whichcommute among themselves, it is an Abelian group as well.The action of the groups G and G on points in space are easily described usingcylindrical coordinates: if x ∈ R is point with cylindrical coordinates ( r, ϑ, z ) , then theaction of the group element Υ mh ∈ G is to send it to a point with cylindrical coordinates ( r, ϑ + 2 πmα, z + mτ ) , while the action of Υ mh • Υ nc ∈ G is to send it to the point withcylindrical coordinates ( r, ϑ + 2 πmα + n Θ , z + mτ ) . In what follows, we will use thenotation (cid:101) Υ to denote a generic isometry from G or G . Given a point x ∈ R , and a group of isometries G (which could be the helical groups G or G described above, for instance), the orbit of x under the group is the set G ◦ x := { (cid:101) Υ ◦ x : (cid:101) Υ ∈ G} . (5)Given a collection of points S ⊂ R and a group of isometries G , we will use the notation G ◦ S to denote the orbits of each of the points in S under the group: G ◦ S := (cid:91) x ∈ S G ◦ x (cid:18) = (cid:91) (cid:101) Υ ∈G (cid:101) Υ ◦ S = (cid:91) (cid:101) Υ ∈G , x ∈ S (cid:101) Υ ◦ x (cid:19) . (6)Let O ⊂ R be a domain with regular boundary that is invariant under a given helicalgroup G , i.e., G ◦ O = O . The symmetry cell or fundamental domain of G in O is a set D ⊂ O such that : (cid:91) (cid:101) Υ ∈G (cid:101) Υ ◦ D = O , (7)and for (cid:101) Υ , (cid:101) Υ ∈ G : ( (cid:101) Υ ◦ D ) (cid:92) ( (cid:101) Υ ◦ D ) = a set of Lebesgue measure for (cid:101) Υ (cid:54) = (cid:101) Υ . (8) rational number. In certain terminology [9, 10], such cases would be identified as rod groups and theterm helical group would be reserved only for cases for which α is an irrational number (i.e., when thegroup is not equivalent to a periodic group generated by a single translation.) However, we will not makethis distinction here. With these hypotheses, the boundary of O , denoted ∂ O , can be shown to be invariant under thegroup as well. In practice, we will require the fundamental domain to have some regularity properties in additionto the conditions in eq. 7 and 8, e.g. it should be connected and have compact closure.
8o see concrete examples of the sets O and D , let D R denote an open disk of radius R on the e , e plane, i.e., D R = { p ∈ R : (cid:107) p (cid:107) R < R } , (9)and let C denote the infinite cylinder obtained by translating D R along e , i.e.: C = D R × { x e : x ∈ R } . (10)Then, the cylinder C has all the properties required of the domain O . Furthermore, weobserve that the finite cylinder D G = D R ×{ x e : 0 ≤ x < τ } serves as the fundamentaldomain of G in C . Finally, the sector described in cylindrical coordinates as: D G = (cid:8) ( r, ϑ, z ) : 0 ≤ r < R, ≤ ϑ < π N , ≤ z < τ (cid:9) , (11)serves as the fundamental domain of G in C . A helical structure i.e., an atomic/molecular structure with helical symmetries issimply the orbit of a set of non-degenerate points under the action of one of the helicalgroups G or G . More precisely, let P G ⊂ D G (cid:0) or P G ⊂ D G in case of G (cid:1) be afinite collection of distinct points labeled (cid:8) x k (cid:9) M G k =1 (cid:0) or (cid:8) x k (cid:9) M G k =1 in case of G (cid:1) . Thesepoints are representative of atomic positions within the fundamental domain and we willrefer to them as simulated points or simulated atoms . The (valence) nuclear chargescorresponding to these atoms will be denoted as (cid:8) Z k (cid:9) M G k =1 (cid:0) or (cid:8) Z k (cid:9) M G k =1 in case of G (cid:1) . Ahelical structure is simply a set of the form: S G , P G = G ◦ P G = (cid:91) (cid:101) Υ ∈G ,k =1 ,...,M G (cid:101) Υ ◦ x k , (12)or S G , P G = G ◦ P G = (cid:91) (cid:101) Υ ∈G ,k =1 ,...,M G (cid:101) Υ ◦ x k . (13)Additionally, for any (cid:101) Υ ∈ G , the atom at the location (cid:101) Υ ◦ x k is taken to be of the samespecies as the atom at x k ∈ P G (similarly also for (cid:101) Υ ∈ G and x k ∈ P G ), and so it isassociated with the same (valence) nuclear charge Z k . The Kohn-Sham equations, as they apply to finite structures can be found in numerousreferences [25, 32, 60]. In order to formulate an appropriate version of the Kohn-Shamequations for helical structures however, we need to keep in mind a few typical featuresof such a structure. In what follows, for the sake of simplicity, we will consider in detailthe case of a structure associated with a helical group generated by a single element(i.e., the group G described above). We will comment on modifications to the abovecase that need to be considered while dealing with a structure associated with a helicalgroup generated two elements (i.e., the group G described above), and present the final9xpressions/equations for this case in Appendix C. A more detailed discussion of themodifications and the application of resulting equations is the scope of future work [70].Helical structures are essentially quasi-one-dimensional in nature. This implies thatthey have limited spatial extent in the e , e plane, while being infinitely extended alongthe e direction. Consequently, it is appropriate to set up the Kohn-Sham equationsfor such a structure in a computational domain which is of limited spatial extent in the e , e plane, while being infinite in extent along e . This, along with the requirementthat a symmetry adapted formulation of the Kohn-Sham equations needs to be solvedon a domain that is also invariant with respect to the symmetry operations of the helicalstructure, suggests the cylinder C as being a natural choice for the computational domain(for a helical structure generated by a single element). The radius of this cylinder hasto be consistent with the requirements that the all the atoms of the helical structureshould be located sufficiently away from the lateral surface of the cylinder so as to allowsufficient decay of various fields that appear in the Kohn-Sham problem.The quasi-one-dimensional nature of the systems under study results in additionalcomplications. Specifically, due to the infinite extent of the system along the e direction,the system is associated with an infinite number of electronic states as well as aninfinite number of nuclei. This potentially poses divergence issues while computing theelectrostatics terms in the Kohn-Sham problem [60, 71] and it is dealt with in this workby solving an appropriate symmetry adapted Poisson problem involving a neutral chargedistribution — such a charge distribution arises as a combination of the electron densityand the nuclear pseudocharges associated with the structure. Additionally, the infinitelymany electronic states have to be incorporated into the Kohn-Sham problem in a mannerthat is consistent with the Pauli exclusion principle and the Aufbau principle [25, 72].Taking cue from the solid state/condensed matter physics literature — specifically, abinitio calculations of crystalline solids [72, 73] — we address this issue here by formulatinga band theory of electronic states for helical structures. This allows the Kohn-Shamproblem for the entire helical structure, as posed on the cylinder C , to be reduced tocomputations on the fundamental domain when augmented with appropriate boundaryconditions.A key ingredient of the band theory for helical structures is an appropriate version ofthe Bloch theorem [37, 73, 74] for such systems. The form of this mathematical result canbe guessed by looking at the analogous case of the Bloch Theorem for one dimensionalperiodic systems and the result appears to have been made use of by earlier researchersin various contexts [8, 50, 63, 64, 66, 75–81]. However, a rigorous mathematical derivationof the result does not seem to appear anywhere in the literature — other than in [8], wherea proof of the existence of Helical Bloch waves was sketched by using tools from the theoryof linear elliptic partial differential equations. In what follows, we address this gap in theliterature and follow up on [8], by establishing the existence and completeness of helicalBloch waves, and then use this to gain insight into the spectrum of the single electronHamiltonian associated with helical systems (i.e., to set up an electronic band theory forsuch systems). This information is subsequently used to set up the governing equationsof the system. Our mathematical treatment closely follows the techniques presented inreferences [8, 82–85]. In general, these states would be expected to be delocalized over the entire volume of the cylinder C . See e.g. equations , in [60]. .2.1. Analysis of the single electron problem for helical structures - helical Bloch waves As a starting point, we consider the single electron Hamiltonian: H = −
12 ∆ + V ( x ) = − (cid:18) ∂ ∂x + ∂ ∂x + ∂ ∂x (cid:19) + V ( x , x , x ) , (14)with the real valued continuous potential V ( x ) invariant under the helical group G , i.e., V ( (cid:101) Υ ◦ x ) = V ( x ) , ∀ (cid:101) Υ ∈ G . (15)This operator naturally arises during each self-consistent field iteration cycle in Kohn-Sham calculations and in that scenario, the invariance of the potential automaticallyfollows from the invariance of the electron density [8].We are interested in functions ψ that satisfy the equation H ψ = λψ within the region C in an appropriate manner. Additionally, to model the decay of the eigenstates asone moves away from the axis of the cylinder to infinity [86, 87], we will enforce Dirichletboundary conditions on the lateral surfaces of the cylinder , i.e., ψ ( x ) = 0 for x ∈ ∂ C . Helical Bloch waves (or helical Bloch states) are solutions to the above equation whichhave the ansatz: ψ ( x ) = e − i πη x τ φ ( x ; η ) , (21) Within the setting of the local density approximation and the use of local pseudopotentials forexample, V ( x ) can be identified as the total effective potential appearing in the Kohn-Sham equationsand can be written as the sum of electrostatic and exchange-correlation terms, i.e., V ( x ) = V es ( ρ ( x )) + V xc ( ρ ( x )) . This “wire” boundary condition is commonly employed in the literature for studying quasi-1D systems[33, 88]. This boundary condition allows the operator H to have some convenient properties withouthaving to enforce any specific decay conditions on V ( x ) as one moves away from the axis of the cylinder. In what follows, we will use the following notation: if A is a measure space with measure µ , then for ≤ p < ∞ , we will use L p ( A, µ ; B ) to denote Lebesgue measurable functions f : A → B which satisfy (cid:90) A (cid:107) f (cid:107) pB dµ < ∞ , and we will use L ∞ ( A, µ ; B ) to denote functions for which ess. sup. x ∈ A (cid:107) f ( x ) (cid:107) B < ∞ .In particular, if A is a domain in R , we will use L ( A ) to denote the usual Hilbert space of complexvalued functions on A which are square integrable (using the Lebsegue measure). Furthermore, H k ( A ) will denote the Sobolev space of tempered distributions whose k th weak derivative lies in L ( A ) , while H ( A ) will denote the subspace of functions in H ( A ) which vanish at the boundary of A in the tracesense. We may view H as an unbounded operator on L ( C ) with the function space Dom. ( H ) = H ( C ) ∩ H ( C ) as the domain of the operator. The operator H is formally symmetric (or, in linear algebra terminology, Hermitian since the underlying function spaces are complex): if f , f are Schwartz functions in C whichobey the boundary condition f ( x ) = f ( x ) = 0 for x ∈ ∂ C , we have: (cid:104) H f , f (cid:105) L ( C ) = − (cid:90) C ∆ f f d x + (cid:90) C V f f d x . (16)On using integration by parts [89] and the decay of f and f as x → ∞ , we get: − (cid:90) C f ∆ f d x = 12 (cid:18) (cid:90) C ∇ f · ∇ f d x − (cid:90) ∂ C f ∇ f · d s (cid:19) . (17)Here d s denotes the oriented surface measure. The second term on the right-hand side above vanishesdue to the boundary conditions obeyed by f , f on ∂ C and so, this leaves us with: (cid:104) H f , f (cid:105) L ( C ) = 12 (cid:90) C ∇ f · ∇ f d x + (cid:90) C V f f d x . (18) φ ( x ; η ) group invariant i.e., φ ( x ; η ) = φ ( (cid:101) Υ ◦ x ; η ) , ∀ (cid:101) Υ ∈ G . (22)and obeys the boundary condition: φ ( x ; η ) = 0 for x ∈ ∂ C , (23)commensurate with the boundary condition on ψ . The parameter η serves a role that isanalogous to k-points in periodic calculations and as shown later, it can be chosen suchthat η ∈ [ − , ) . In what follows, we first show the existence of such solutions and thendemonstrate their completeness. In essence, these results together give us informationthat certain special electronic states (i.e., helical Bloch states) can be always found to beassociated with the single electron Hamiltonian of a helical structure, and they furtherinform us that such special states can be used to characterize all of the possible electronicstates of the system (within the single electron model). Therefore, it is sufficient for usto restrict our attention to these states.First, to demonstrate the existence of these special solutions, we have: Theorem 2.3 ( Existence theorem for helical Bloch waves ) . Let V ( x ) be a realvalued continuous potential that is invariant under the helical group G and let H denotethe operator − ∆ + V ( x ) . For any η ∈ R there exist a countable number of solutionsof the equation H ψ = λψ in C which are expressible in terms of the helical Bloch ansatz(eqs. 21-23).Proof. We fix η ∈ R and substitute the helical Bloch wave ansatz in the equation H ψ = λψ to find that φ ( x ; η ) should obey the following auxiliary equation in the region C : h aux η φ = − (cid:18) ∆ φ − i πητ ∂φ∂x − π η τ φ (cid:19) + V φ = λφ . (24)Additionally, φ should be group invariant and obey the zero Dirichlet boundary conditionon ∂ C . Let D denote the interior of the fundamental domain D G , i.e., it is the open setdescribed in cylindrical coordinates as D = { ( r, ϑ, z ) : 0 ≤ r < R, < z < τ } . Theboundary of D includes the lateral surface ∂ D r = R = { ( r, ϑ, z ) : r = R, < z < τ } that is shared with C , as well as the discs ∂ D z =0 = { ( r, ϑ, z ) : 0 ≤ r ≤ R, z = 0 } and ∂ D z = τ = { ( r, ϑ, z ) : 0 ≤ r ≤ R, z = τ } , which are both parallel to the e − e plane.The group operation Υ h (i.e., the generator of the group G ) maps ∂ D z =0 to ∂ D z = τ andconversely, Υ − h maps ∂ D z = τ to ∂ D z =0 .We now restrict the auxiliary eigenvalue problem h aux η φ = λφ as outlined in eq. 24,to the region D by imposing the boundary conditions φ ( x , η ) = φ (Υ h ◦ x , η ) , ∇ φ ( x , η ) = In a similar manner, we get: (cid:104) f , H f (cid:105) L ( C ) = 12 (cid:90) C ∇ f · ∇ f d x + (cid:90) C f V f d x (19) = (cid:104) H f , f (cid:105) L ( C ) , (20)as the potential V ( x ) is real. Since Schwartz functions are dense in the domain of H , the result follows.The direct integral decomposition of H (Appendix B) makes it easy to appreciate that H is in factself-adjoint. − πα ∇ φ (Υ h ◦ x , η ) for x ∈ ∂ D z =0 and (as before), φ ( x , η ) = 0 for x ∈ ∂ D r = R . Theoperator h aux η on L ( D ) is uniformly elliptic, and as shown in Appendix A, it is alsosymmetric with the above boundary conditions. Since D is a bounded domain and V ( x ) ∈ L ∞ ( D ) , the operator h aux η satisfies the conditions of Gårding’s Inequality (Theorems 9.17,9.18 in [90]; Section 6.2 in [89]). This guarantees that h aux η has a unique self-adjointextension in L ( D ) , which we also denote as h aux η here. Furthermore, as a consequenceof the Rellich-Kondrachov Compactness Theorem (Theorem 7.29 in [90]; Section 5.7 in[89]), h aux η can be shown to have a compact resolvent (Lemma 9.20 in [90]). Consequently, h aux η has a discrete set of eigenvalues λ j ( η ) and corresponding eigenfunctions φ j ( x ; η ) (Theorem 6.29 in [91]; Theorem 9.22 in [90]). Each eigenvalue is of finite multiplicity andsuch that λ j ( η ) → ∞ as j → ∞ . Results from elliptic regularity theory (Sections 9.5,9.6 in [90]; Section 6.3 in [89]) imply that φ j ( · ; η ) ∈ H ( D ) . We now use the boundaryconditions on φ outlined above to extend the eigenfunctions φ j ( x ; η ) to all of C , notingthat these boundary conditions are meaningful in the trace sense since the eigenfunctionsare in H ( D ) . Thereafter, defining ψ j ( x , η ) = e − i πη x τ φ j ( x , η ) , for j ∈ N and x ∈ C ,establishes the theorem. (cid:4) We define
Λ = { λ j ( η ) : η ∈ R , j ∈ N } and Ψ = { ψ j ( · ; η ) : η ∈ R , j ∈ N } as thecollection of generalized eigenvalues and generalized eigenfunctions associated with H .The first observation we make is that the sets Λ and Ψ are unchanged upon restricting η ∈ [ − , ) . To see this, we recall that λ j ( η ) and ψ j ( · ; η ) are obtained by computingthe spectrum of H when subjected to the conditions in eqs. 21, 22. However, theseequations can be equivalently recast as the following condition on ψ : ψ (Υ h ◦ x ) = e − i πη x ττ φ (Υ h ◦ x ; η ) = e − i πη e − i πη x τ φ ( x ; η ) = e − i πη ψ ( x ) (25)or more generally, for m ∈ Z : ψ (Υ m h ◦ x ) = e − i πmη ψ ( x ) . (26)In other words, solving H ψ = λψ while imposing the condition ψ (Υ h ◦ x ) = e − i πη ψ ( x ) on ψ also gives us the sets Λ and Ψ . Since e − i πη = e − i π ( η + n ) for any n ∈ Z , we see thatthe boundary conditions on ψ do not change upon translating the value η by an integer.Thus, it suffices to restrict η ∈ [ − , ) . In what follows, we will denote I = [ − , ) , andwe will re-define Λ = { λ j ( η ) : η ∈ I , j ∈ N } and Ψ = { ψ j ( · ; η ) : η ∈ I , j ∈ N } . In keepingwith solid state physics terminology, we will refer to the set I as reciprocal space (or morespecifically, the Brillouin zone of the reciprocal space). Consequently, the dependence ofa quantity on η will be termed as reciprocal space dependence while its dependence onusual physical space will be termed as real space dependence.For a given j ∈ N , we will refer to the set Λ j = { λ j ( η ) : η ∈ I ) } as a helical band .Results from the theory of regular perturbations of self-adjoint problems [91, 92] imply The real numbers λ j ( η ) are generalized eigenvalues of H since (as discussed later) they are part ofthe essential spectrum of H and not its discrete spectrum. On a similar note, the functions ψ j ( · ; η ) do notbelong in L ( C ) and therefore, they are not eigenfunctions of H in the usual sense. However, as discussedabove, they do satisfy an equation of the form H ψ j ( · ; η ) = λ ( η ) ψ j ( · ; η ) , thus suggesting their similarityto conventional eigenvalues and eigenfunctions. The Dirichlet boundary condition in eq. 23 is also obeyed equivalently by ψ and does not need to befurther considered here. j ∈ N ) the map η (cid:55)→ λ j ( η ) is analytic. Therefore, the set Λ j is connectedand compact. We will refer to the set
Ψ = { ψ j ( · ; η ) : η ∈ I , j ∈ N } as the collection of helical Blochstates corresponding to the helical bands. If we fix η ∈ I , then the set Ψ η = { ψ j ( · ; η ) : j ∈ N } has the property that it is orthonormal and complete in L ( D ) . This follows directlyfrom the properties of the group invariant functions φ j ( x ; η ) defined above. Specifically,for j, j (cid:48) ∈ N : (cid:104) ψ j ( · ; η ) , ψ j (cid:48) ( · ; η ) (cid:105) L ( D ) = (cid:90) D ψ j ( x , η ) ψ j (cid:48) ( x , η ) d x = (cid:90) D e − i πη x τ φ j ( x , η ) e i πη x τ φ j (cid:48) ( x , η ) d x = (cid:90) D φ j ( x , η ) φ j (cid:48) ( x , η ) d x = δ j,j (cid:48) . (28)Furthermore, if h ∈ L ( D ) such that (cid:104) h ( · ) , ψ j (cid:48) ( · ; η ) (cid:105) L ( D ) = 0 for every j ∈ N , then wemust have: (cid:90) D h ( x ) ψ j ( x , η ) d x = (cid:90) D h ( x ) e i πη x τ φ j ( x , η ) d x = 0 , ∀ j ∈ N . (29)Due to the completeness of the functions φ j ( · ; η ) it then follows that h ( x ) e i πη x τ = 0 ,i.e., h ( x ) = 0 almost everywhere in D . Thus, the set Ψ η is complete in L ( D ) .Due to the completeness of the set Ψ η for each η ∈ I , it actually follows that theset of helical Bloch states (i.e., the set Ψ ) is complete in L ( C ) . To prove this importantresult, we first need to establish a few preliminaries related to the so-called Bloch-Floquettransform [83, 85], as we now do. Lemma 2.4.
Let f ∈ L ( C ) , η ∈ I and m ∈ Z . We define: g ( x , η ) := (cid:88) m ∈ Z f (Υ m h ◦ x ) e i πmη . (30) Then g is defined almost everywhere in D and further, g ∈ L ( D × I ) . In contrast to the rigorous proof presented above, a formal derivation of the Bloch theorem for ahelical structure, inspired by the solid state physics literature [37, 38] is as follows: We observe thatsince the Laplacian commutes with all isometry operations – including those that constitute the group G , and further, since the potential V ( x ) is group invariant (eq. 15), the operator H must commute withthe symmetry operations in the group G . Specifically, for any continuous function f defined over C , wemay define the operators: T = (cid:8) T (cid:101) Υ : T (cid:101) Υ f ( x ) = f ( (cid:101) Υ − ◦ x ) (cid:9) (cid:101) Υ ∈G . (27)Then, for any function f in the domain of H , the relationship T (cid:101) Υ H f = H T (cid:101) Υ f holds for any T (cid:101) Υ ∈ T . Thiscommutation property can be used to infer that the unitary representations of G and the operator H can be “simultaneously diagonalized” in a suitable basis of common “eigenstates”. Since G is an Abeliangroup, its irreducible representations are all one-dimensional [41, 93]. Furthermore, these irreduciblerepresentations can be used to decompose any unitary representation of the group [40, 41]. This suggeststherefore that the eigenstates associated with H transform under the group in a manner similar tothe irreducible representations of G , which then implies the helical Bloch theorem. While the aboveargument is perhaps correct in spirit and variants of the argument appear often in the physics literature(in the context of periodic systems) it has a number of technical deficiencies owing to the fact that H isan unbounded operator and the group G is infinite. These issues prevent heuristic arguments like theone above – which are more suited to representations of finite groups on finite dimensional spaces – frombeing applied in the current context. In particular e.g., Bloch states are not eigenfunctions in the usualsense since they are not square integrable. roof. We denote ˜ f m ( x ) = f (Υ m h ◦ x ) . Then, g ( x , η ) = (cid:88) m ∈ Z ˜ f m ( x ) e i πmη . To demonstratethat g is defined almost everywhere in D , It suffices to show that for almost every x ∈ D ,the sequence (cid:8) ˜ f m ( x ) (cid:9) m ∈ Z is square summable, since in that case the expression for g canbe simply interpreted as a Fourier series expansion [94]. To this end, we have: (cid:88) m ∈ Z | ˜ f m ( x ) | ≤ (cid:90) x ∈D (cid:88) m ∈ Z | ˜ f m ( x ) | d x = (cid:88) m ∈ Z (cid:90) x ∈D | ˜ f m ( x ) | d x = (cid:88) m ∈ Z (cid:90) x ∈D | f (Υ m h ◦ x ) | d x = (cid:88) m ∈ Z (cid:90) x ∈ Υ m h ◦D | f ( x ) | d x = (cid:90) x ∈C | f ( x ) | d x < ∞ , (31)since f ∈ L ( C ) . The exchange of the integral and the summation is justified using theFubini-Tonelli Theorem [94]. Since this establishes (cid:8) ˜ f m ( x ) (cid:9) m ∈ Z as a square summablesequence, we may now use Parseval’s identity [94] and eq. 30 to obtain: (cid:90) I | g ( x , η ) | dη = (cid:88) m ∈ Z | ˜ f m ( x ) | . (32)Integrating both sides of this expression for x ∈ D and using the steps in eq. 31 establishesthat g ∈ L ( D × I ) . (cid:4) The following result is the converse of Lemma 2.4 and is established using the sametools as above:
Lemma 2.5.
Let g ∈ L ( D × I ) and for x ∈ D , m ∈ Z , let: ˜ f m ( x ) = (cid:90) I g ( x , η ) e − i πmη dη . (33) Furthermore, let the function f be an extension of ˜ f from the domain D to the domain C in the sense that for x ∈ D , f (Υ m h ◦ x ) := ˜ f m ( x ) . (34) Then, f ∈ L ( C ) .Proof. By Tonelli’s theorem [94], since g ∈ L ( D × I ) , it holds that g ( x , · ) ∈ L ( I ) foralmost every x ∈ D . Then, we may interpret eq. 33 as a Fourier transform in η . ByParseval’s identity [94], we have: (cid:90) I | g ( x , η ) | dη = (cid:88) m ∈ Z | ˜ f m ( x ) | . (35)Integrating both sides over x ∈ D and using g ∈ L ( D × I ) , we get: ∞ > (cid:90) x ∈D (cid:18) (cid:90) I | g ( x , η ) | dη (cid:19) d x = (cid:90) x ∈D (cid:18) (cid:88) m ∈ Z | ˜ f m ( x ) | (cid:19) d x = (cid:88) m ∈ Z (cid:90) x ∈D (cid:88) m ∈ Z | ˜ f m ( x ) | d x = (cid:88) m ∈ Z (cid:90) x ∈D | f (Υ m h ◦ x ) | d x = (cid:88) m ∈ Z (cid:90) x ∈ Υ m h ◦D | f ( x ) | d x = (cid:90) x ∈C | f ( x ) | d x . (36)15his shows that f ∈ L ( C ) , as required. Note that the interchange of the summation andthe integral in the calculations above can be justified using the Fubini-Tonelli Theorem[94]. (cid:4) Lemma 2.4 establishes the existence of an operator U : L ( C ) → L ( D × I ) defined as: ( U f )( x , η ) = (cid:88) m ∈ Z f (Υ m h ◦ x ) e i πmη , (37)while Lemma 2.5 establishes the existence of its inverse U − : L ( D × I ) → L ( C ) definedas: ( U − g )(Υ m h ◦ x ) = (cid:90) I g ( x , η ) e − i πmη dη . (38)To verify that eq. 38 indeed defines the inverse of the operator in eq. 37, we consider f ∈ L ( C ) and g ∈ L ( D × I ) such that g = U f , i.e., g ( x , η ) = (cid:88) m ∈ Z f (Υ m h ◦ x ) e i πmη . (39)We now multiply the above by e − i πm (cid:48) η for m (cid:48) ∈ Z and integrate over η , to arrive at: (cid:90) I g ( x , η ) e − i πm (cid:48) η dη = (cid:90) I (cid:88) m ∈ Z f (Υ m h ◦ x ) e i π ( m − m (cid:48) ) η dη = (cid:88) m ∈ Z (cid:90) I f (Υ m h ◦ x ) e i π ( m − m (cid:48) ) η dη = f (Υ m (cid:48) h ◦ x ) . (40)Thus f = U − g in accordance with eq. 38.We also observe, based on the calculations in eq. 36 that: (cid:107) f (cid:107) L ( C ) = (cid:107)U f (cid:107) L ( D× I ) , (41)and therefore, the operator U is an isometric-isomorphism between the spaces L ( C ) and L ( D × I ) . In analogy to the Bloch-Floquet transform in the literature used for studyingperiodic problems [83, 85], we will refer to the operator U as the helical Bloch-Floquettransform . This operator allows us to demonstrate the completeness of the helicalBloch waves in L ( C ) . The basic idea behind this proof is to map a given f ∈ L ( C ) to itscounterpart in L ( D × I ) and to then use the completeness of the set Ψ η for each η ∈ I . Theorem 2.6 ( Completeness theorem for helical Bloch waves ) . Let f ∈ L ( C ) ,and for (cid:96) ∈ N , x ∈ C , let: f (cid:96) ( x ) := (cid:96) (cid:88) s =1 (cid:90) I (cid:10) ( U f )( · ; η ) , ψ s ( · ; η ) (cid:11) L ( D ) ψ s ( x , η ) dη . (42) Then f (cid:96) → f in L ( C ) as (cid:96) → ∞ . Eq. 41 shows that U is an isometry and Lemma 2.5 shows that it has a well defined inverse. Therefore,it is a unitary operator [94]. This operator is closely related to the so-called Zak transform [95, 96] associated with the group. roof. Since U f ∈ L ( D × I ) , it follows from Fubini’s theorem that U f ( · ; η ) ∈ L ( D ) foralmost every η ∈ I . Therefore, it can be approximated using the functions in the set Ψ η .Consequently, if we define: g (cid:96) ( x , η ) := (cid:96) (cid:88) s =1 (cid:10) ( U f )( · ; η ) , ψ s ( · ; η ) (cid:11) L ( D ) ψ s ( x , η ) , (43)then g (cid:96) ( · ; η ) → ( U f )( · ; η ) in L ( D ) as (cid:96) → ∞ for almost every η ∈ I . In other words, theresidual: r (cid:96) ( η ) = (cid:107) ( U f )( · ; η ) − g (cid:96) ( · ; η ) (cid:107) L ( D ) , (44)has the property that r (cid:96) ( η ) → for almost every η ∈ I , as (cid:96) → ∞ . Furthermore usingthe identity (cid:107) f + f (cid:107) ≤ (cid:107) f (cid:107) + (cid:107) f (cid:107) ) , as well as Bessel’s inequality [94] on g (cid:96) → U f ,we get: r (cid:96) ( η ) = (cid:107) ( U f )( · ; η ) − g (cid:96) ( · ; η ) (cid:107) L ( D ) ≤ (cid:18) (cid:107) ( U f )( · ; η ) (cid:107) L ( D ) + (cid:107) g (cid:96) ( · ; η ) (cid:107) L ( D ) (cid:19) ≤ (cid:18) (cid:107) ( U f )( · ; η ) (cid:107) L ( D ) + (cid:107) ( U f )( · ; η ) (cid:107) L ( D ) (cid:19) = 4 (cid:107) ( U f )( · ; η ) (cid:107) L ( D ) . (45)However, (cid:107) ( U f )( · ; η ) (cid:107) L ( D ) is in L ( I ) based on the calculations in Lemma 30. Therefore,by the Dominated Convergence Theorem [94]: (cid:90) I r (cid:96) ( η ) dη → as (cid:96) → ∞ , (46)and consequently: (cid:107)U f − g (cid:96) (cid:107) L ( D× I ) → as (cid:96) → ∞ . (47)Since U is an isometric isomorphism, this implies that U − g (cid:96) → f in L ( C ) as (cid:96) → ∞ .Now, using eq. 38, we see that: U − g (cid:96) (Υ m h ◦ x ) = (cid:90) I (cid:96) (cid:88) s =1 (cid:10) ( U f )( · ; η ) , ψ s ( · ; η ) (cid:11) L ( D ) ψ s ( x , η ) e − i πmη dη . (48)On the other hand, evaluating eq. 42 at x = Υ m h ◦ y , and using eq. 26, we see that: f (cid:96) (Υ m h ◦ y ) = (cid:96) (cid:88) s =1 (cid:90) I (cid:10) ( U f )( · ; η ) , ψ s ( · ; η ) (cid:11) L ( D ) ψ s (Υ m h ◦ y , η ) dη (49) = (cid:96) (cid:88) s =1 (cid:90) I (cid:10) ( U f )( · ; η ) , ψ s ( · ; η ) (cid:11) L ( D ) e − i πmη ψ s ( y , η ) dη . (50)Comparing eqs. 48 and 50, it follows that U − g (cid:96) = f (cid:96) since y and m are generic, andtherefore, f (cid:96) → f in L ( C ) when (cid:96) → ∞ , as required. (cid:4)
17s mentioned earlier, the above results imply in essence that the spectral properties of H can be described completely in terms of helical bands and helical Bloch states (refer toAppendix B for further discussion along these lines). Additionally, since the behaviorof any helical Bloch state over all of C is completely specified based on its behaviorover D (once a value of η ∈ I is chosen), the single electron problem posed on all of C can be reduced to a set of problems (indexed by η ) posed on the fundamental domain(illustrated in Figure 2). Consequently, by appropriate use of the helical Bloch states andthe helical bands, quantities of interest in Kohn-Sham theory (which can be describedusing the solutions to the single electron problem), can be formulated entirely in terms ofquantities specified on the fundamental domain. We now look at this procedure in moredetail. In what follows, we will consider the helical structure to be at finite electronic tem-perature T e Kelvin and we will ignore spin polarization effects. For the sake of clarity ofpresentation, we will itemize the formulation/derivation of the various terms and equa-tions, as we go along.Electron Density and Density Matrix: A quantity of key importance in Kohn-Shamtheory is the electron density ρ ( x ) . For a finite structure, such as a molecule or a cluster, An immediate consequence of the completeness theorem for Bloch states is that the spectrum of H is completely contained in the set of helical bands, i.e., more precisely, spec . ( H ) ⊆ clos . (Λ) , with clos . ( · ) denoting the (topological) closure. This is because, if κ ∈ R is such that κ / ∈ clos . (Λ) , then the action of ( H − κ ) − on f ∈ L ( C ) can be computed formally using eq. 42 in Theorem 2.6 as: ( H − κ ) − f = ∞ (cid:88) s =1 (cid:90) I (cid:10) ( U f )( · ; η ) , ψ s ( · ; η ) (cid:11) L ( D ) ( H − κ ) − ψ s ( · ; η ) dη . (51)Now, using H ψ s ( · ; η ) = λ s ( η ) ψ s ( · ; η ) , we have: ( H − κ ) − ψ s ( · ; η ) = 1 λ s ( η ) − κ ψ s ( · ; η ) , (52)so that: ( H − κ ) − f = ∞ (cid:88) s =1 (cid:90) I (cid:10) ( U f )( · ; η ) , ψ s ( · ; η ) (cid:11) L ( D ) λ s ( η ) − κ ψ s ( · ; η ) dη . (53)Since κ / ∈ clos . (Λ) , the term λ s ( η ) − κ remains bounded even as s → ∞ . Therefore, the right-hand sideof eq. 53 can be interpreted as a bounded operator on f , and so, κ must belong to the resolvent set of H . Conversely, based on the techniques presented in [82, 83] it is also possible to directly demonstratethat Λ ⊆ spec . ( H ) , by constructing a suitable singular sequence of the form u l ( x ) = ζ (cid:18) x l (cid:19) ψ j ( x , η ) ,and using Weyl’s criterion [91, 92, 97]. Here ζ ( · ) is a carefully chosen smooth cutoff function. Since, spec . ( H ) is always a closed set [91], and by definition, clos . (Λ) is the smallest closed set containing Λ , itfollows that clos . (Λ) = spec . ( H ) . Furthermore, if λ ∈ spec . ( H ) , it can be immediately seen to be partof the essential spectrum of H . This is because if it were part of the point spectrum, then λ would beassociated with eigenfunctions of finite multiplicity. Due to the fact that H commutes with the operatorsin T = (cid:8) T (cid:101) Υ : T (cid:101) Υ f ( x ) = f ( (cid:101) Υ − ◦ x ) (cid:9) (cid:101) Υ ∈G , these eigenfunctions would be left invariant by the operators in T as well (also see footnote 24). However, this would contradict the requirement that these eigenfunctionsbelong to L ( C ) . Note that these above results also follow from the direct integral decomposition of theHamiltonian discussed in Section 2.2.2 and Appendix B. elical Bloch TheoremSingle electron problem posed on infinite cylinder Single electron problem reduced to fundamental domain 𝒞 e e Dirichlet boundary conditionsDirichlet boundary conditions Helical Bloch boundary conditionsHelical Bloch boundary conditions ∂𝒟 z =0 ∂𝒟 z = τ ∂𝒟 r = R ∂𝒞 𝒟 Helical structure
Figure 2: Illustration of the Helical Bloch Theorem (Theorem 2.3) for the case when the single electronHamiltonian obeys the symmetries of the helical group G . this can be expressed in a straightforward manner in terms of the associated Kohn-Shameigenstates and electronic occupations [32, 60]. For a helical structure however, care hasto be taken to express this quantity due to the fact that there are effectively an infinitenumber of electrons associated with the structure. In what follows, motivated by rigorousmathematical results related to the description of electronic states in crystalline systems[98–100], we address this issue by first defining the single particle density operator [101–103] in terms of the single electron Hamiltonian: ˜Γ = f T e ( H ) , (54)and then expressing the electron density in terms of the diagonal of this operator, i.e., incoordinate representation: ρ ( x ) = 2 ˜Γ( x , x ) . (55)In eq. 54 above, f T e ( · ) denotes the Fermi-Dirac distribution function at electronic tem-perature T e : f T e ( y ) = 11 + exp (cid:0) y − λ F k B T e (cid:1) , (56)and λ F and k B denote the Fermi level and the Boltzmann constant respectively.In order to develop a more transparent expression for the electron density in accor-dance with eq. 55, it is useful to first recast eq. 54 in terms of helical Bloch states. Asshown in Appendix B, the apparatus of direct integrals [85, 104] allows us to do this in19 mathematically rigorous manner. The key result from Appendix B is that the helicalBloch-Floquet transform allows the single electron Hamiltonian to be “block-diagonalized”into a set of problems associated with the helical Bloch states that are posed over thefundamental domain, i.e.: U H U − = (cid:90) ⊕ I H η dη . (57)Here, as before, H represents the operator − ∆ + V ( x ) over the cylinder C along withthe boundary condition ψ ( x ) = 0 for x ∈ ∂ C . The potential V ( x ) is group invariant(eq. 15) and the unitary operator U : L ( C ) → L ( D × I ) represents the helical Bloch-Floquet transform (eq. 37). The operators { H η } η ∈ I represent the fibers of H (in the senseof direct integrals) and are closely related to the operators h aux η introduced in the proofof Theorem 2.3 (eq. 24) . Specifically, for each η ∈ I , the operator H η represents theoperator − ∆ + V ( x ) over the interior of the fundamental domain (i.e., the set D ) alongwith the boundary conditions ψ (Υ h ◦ x ) = e − i πη ψ ( x ) , R − πα ∇ ψ (Υ h ◦ x ) = e − i πη ∇ ψ ( x ) for x ∈ ∂ D z =0 and, ψ ( x ) = 0 for x ∈ ∂ D r = R . The eigenstates of the operators { H η } η ∈ I are precisely the helical bands Λ = { λ j ( η ) : η ∈ I , j ∈ N } and the helical Bloch states Ψ = { ψ j ( · ; η ) : η ∈ I , j ∈ N } restricted to the region D . We may therefore, express theoperators H η using a spectral representation [91, 105, 106] as: H η = ∞ (cid:88) j =1 λ j ( η ) ψ j ( · ; η ) ψ j ( · ; η ) (58)Using the above expression, noting that the function f T e ( · ) is analytic [107] and finally,also making use of the properties of the direct integral representation [85, TheoremXIII.85], we get that as an operator on L ( C ) , the density matrix admits the repre-sentation: ˜Γ = U − (cid:18) (cid:90) ⊕ I f T e ( H η ) dη (cid:19) U = U − (cid:18) (cid:90) ⊕ I f T e (cid:0) λ j ( η ) (cid:1) ψ j ( · ; η ) ψ j ( · ; η ) dη (cid:19) U . (59)In other words, as an operator on L ( D × I ) (which is more appropriate and convenient todeal with, while describing quantities over the fundamental domain), the density matrixadmits the following expression in coordinate representation (with x , y ∈ D ): Γ( x , y ) = (cid:90) ⊕ I ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1) ψ j ( x ; η ) ψ j ( y ; η ) dη . (60)From this, it follows that the electron density can be expressed as (for x ∈ D ): ρ ( x ) = 2 (cid:90) I ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1) | ψ j ( x ; η ) | dη . (61) The key difference is that the operators h aux η include η dependence in the operators themselves andhave group invariant solutions, whereas the operators H η include η dependence in the boundary conditionsand have helical Bloch solutions (i.e., solutions which are group invariant up to an η dependent phase).
20t is easy to see from the above expression that the electron density is group invariantand also obeys a zero-Dirichlet boundary condition on the lateral surface of D (i.e., for x ∈ ∂ D r = R ). It is also apparent from the expression for the density matrix (eq. 60) thatthe following invariance relationship holds for any m ∈ Z and x , y ∈ D : Γ( x , y ) = Γ(Υ m h ◦ x , Υ m h ◦ y ) . (62)For notational simplicity, it is convenient to introduce the scalars g j ( η ) = f T e (cid:0) λ j ( η ) (cid:1) ,i.e. the thermalized occupation numbers of the electronic states of the system. We willdenote the collection of occupation numbers as G = (cid:8) g j ( η ) = f T e (cid:0) λ j ( η ) (cid:1) : η ∈ I , j ∈ N } .The electron density for an extended system is expected to obey the constraint ofhaving a fixed number of electrons per unit fundamental domain of the system, eventhough the electronic states themselves are delocalized over the entire structure [25, 43].Denoting the number of electrons per unit cell as N e , in our case, this leads to: (cid:90) D ρ ( x ) d x = N e , (63)from which, using the orthonormality of the Bloch states over D , follows the constraint: (cid:90) I ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1) dη = 2 (cid:90) I ∞ (cid:88) j =1 g j ( η ) = N e . (64)In practice, the above equation can be used to compute the Fermi-level ( λ F ) of the system.In what follows, we will sometimes need to express eq. 59 in a manner that is similarto eq. 57, i.e., Γ = U ˜Γ U − = (cid:90) ⊕ I ˜Γ η dη , (65)with the fibers of ˜Γ denoted as: ˜Γ η = ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1) ψ j ( x ; η ) ψ j ( y ; η ) = ∞ (cid:88) j =1 g j ( η ) ψ j ( x ; η ) ψ j ( y ; η ) . (66)Electronic Free Energy: With the above expressions in place, we can now use anenergy-minimization formalism to deduce the governing equations of Kohn-Sham theoryfor the helical structure. Since the structure is infinite, the quantity of primary impor-tance in this regard is the electronic free energy per unit fundamental domain , denotedhere as F (Λ , Ψ , P G , D , G ) . This notation emphasizes the dependence of this quantityon the helical Bloch states, the helical Bloch bands, the positions of the simulated atoms P G = (cid:8) x k (cid:9) M G k =1 within the fundamental domain, (the interior of) the fundamental do-main and the helical group itself. Following [60], we express this quantity within thepseudopotential [108, 109] and Local Density Approximations [24] as: F (Λ , Ψ , P G , D , G ) = T s (Λ , Ψ , D , G ) + E xc ( ρ, D ) + K (Λ , Ψ , P G , D , G )+ E el ( ρ, P G , D , G ) − T e S (Λ) , (67) We would like to thank Eric Cances (Ecole des Ponts ParisTech) and Carlos Garcia Cervera (Univ.of California, Santa Barbara) for email communication related to technicalities of the above derivationof eq. 61 and also for providing useful references. R ) with single particle density matrix D ,the kinetic energy can be expressed as [98, 102, 103]: T finites = 2 Tr. (cid:2) −
12 ∆ D (cid:3) , (68)with Tr. [ · ] denoting the operator trace (of a trace-class operator on L ( R ) ). Analogously,it would make sense to consider the trace per unit fundamental domain in case of thehelical structure. As described in Appendix B, for an operator A which is invariantunder the group G and which is locally trace-class, it is possible to assign meaning tothe trace per unit fundamental domain by means of the direct integral decomposition.Specifically, if the helical Bloch-Floquet transform block-diagonalizes the operator intoits fibers as: U A U − = (cid:90) ⊕ I A η dη , (69)then the trace per unit cell (denoted Tr. [ · ] henceforth) can be expressed as:Tr. [ A ] = (cid:90) I Tr. [ A η ] dη , (70)with the trace inside the integral signifying the usual operator trace in L ( D ) . Theexpression for the kinetic energy per unit fundamental domain for the helical structuretherefore boils down to: T s (Λ , Ψ , D , G ) = 2 Tr. (cid:2) −
12 ∆ ˜Γ (cid:3) = 2 (cid:90) I Tr. (cid:20)(cid:0) −
12 ∆ ˜Γ (cid:1) η (cid:21) dη . (72)We now write (cid:0) − ∆˜Γ (cid:1) η as ( − ∆) η ˜Γ η , observe that ˜Γ η is already available from eq. 66.Next, based on the discussion in Appendix B, we note that the fibers of the Laplacianon L ( C ) are simply the Laplacian operators on L ( D ) with ( η -dependent) helical Bloch The operator trace for any trace-class operator O on L ( D ) can be computed as:Tr. [ O ] = ∞ (cid:88) j =1 (cid:104) O f j , f j (cid:105) L ( D ) , (71)where { f j } ∞ j =1 can be any orthonormal basis of L ( D ) . Refer e.g. to [98, 102] for broader discussions oftrace-class and locally trace-class operators in the context of electronic structure models. Since the helical Bloch states ψ j ( · ; η ) belong to the domain of the operator H η , it follows that thetraces in eq. 72 are finite, making the expressions in that equation well defined. See [98, 102] for furthermathematical details along these lines. ψ j ( x ; η ) already satisfy these boundary conditions,and they form a basis of L ( D ) , it follows that: T s (Λ , Ψ , D , G ) = 2 (cid:90) I (cid:20) ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1) (cid:90) D −
12 ∆ ψ j ( x ; η ) ψ j ( y ; η ) d x (cid:21) dη = (cid:90) I ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1)(cid:10) ∆ ψ j ( · ; η ) , ψ j ( · ; η ) (cid:11) L ( D ) dη = (cid:90) I ∞ (cid:88) j =1 g j ( η ) (cid:10) ∆ ψ j ( · ; η ) , ψ j ( · ; η ) (cid:11) L ( D ) dη . (73)Exchange-Correlation Term: The second term on the right-hand side of eq. 67 rep-resents the exchange correlation energy of the electrons per unit fundamental domain.Within the Local Density Approximation (LDA) [24], it can be written as: E xc ( ρ, D ) = (cid:90) D ε xc [ ρ ( x )] ρ ( x ) d x . (74)Note that it is also possible to modify this expression to use more sophisticated exchangecorrelation functionals such as the Generalized Gradient Approximation [110] and thiswill have little bearing on our subsequent discussion.Nonlocal Pseudopotential Energy Term: The third term on the right-hand side ofeq. 67 represents the energetic contribution from the nonlocal part of the pseudopo-tential. For a finite system of N at atoms located at the points { p k ∈ R } N at k =1 this termhas the following form: K finite = 2 Tr. (cid:2) V nl D (cid:3) , (75)with the non-local pseudoptential operator V nl expressible in Kleinman-Bylander form[111] as: V nl = N at (cid:88) k =1 (cid:88) p ∈N k γ k,p χ k,p ( · ; p k ) χ k,p ( · ; p k ) . (76)Here, N k denotes the collection of projectors associated with the atom at p k , χ k,p arethe projection functions, and γ k,p are the corresponding normalization constants. Thefunctions χ k,p are themselves expressible in terms of atomic orbitals and are usuallysupported in a small region of space by design [112]. To obtain the correct analog ofthis expression for the helical structure i.e., the nonlocal pseudopotential energy per unitfundamental domain , it is useful to recall that this contribution to the energy is tied tothe atoms in the fundamental domain as well as the electronic states in the system. Itis of a somewhat different nature as compared to the kinetic energy term for instance,in which case the electrons are the only contributing source. Since the electrons in theextended structure are delocalized, the trace per unit fundamental domain leads to theappropriate expression in that case. In case of the nonlocal pseudopotential energy termhowever, the contribution from the electrons is delocalized, while those from the atoms23re not. With this in mind, we now focus on the atoms in the fundamental domain,and denote the non-local pseudoptential operator associated with these atoms as: V nl D = M G (cid:88) k =1 (cid:88) p ∈N k γ k,p χ k,p ( · ; x k ) χ k,p ( · ; x k ) . (77)Then, in analogy with eq. 75, the nonlocal pseudopotential energy per unit fundamentaldomain in case of the helical structure can be written by considering the action of V nl D tothe density matrix operator defined in eq. 59, i.e.: K = 2 Tr. [ V nl D ˜Γ ] . (78)To simplify this expression , we use eq. 65, employ the unitarity of the operator U , aswell as the invariace of the trace under unitary transformations to obtain: K = 2 Tr. [ V nl D U − Γ U ] = 2 Tr. [ U − U V nl D U − Γ U ] = 2 Tr. [ U V nl D U − Γ ] . (79)Next, we observe that: U V nl D U − = U (cid:18) M G (cid:88) k =1 (cid:88) p ∈N k γ k,p χ k,p ( · ; x k ) χ k,p ( · ; x k ) (cid:19) U − = M G (cid:88) k =1 (cid:88) p ∈N k γ k,p U χ k,p ( · ; η ; x k ) U χ k,p ( · ; η ; x k ) . (80)In what follows, for the sake of brevity, we will denote the helical Bloch-Floquet transformof the projection functions, i.e., U χ k,p ( · ; η ; x k ) as ˆ χ k,p ( · ; η ; x k ) and note that they can berepresented via eq. 37 (for x ∈ D and η ∈ I as): ˆ χ k,p ( x ; η ; x k ) = U χ k,p ( x ; η ; x k ) = (cid:88) m ∈ Z χ k ; p (Υ m h ◦ x ; x k ) e i πmη . (81)In practice, since the projection functions often have small support (centered about atomicpositions), it is possible to truncate the above summation to just a few terms. With We would like to thank Phanish Suryanarayana, Georgia Institute of Technology, for discussionswhich helped clarify some of the properties of the nonlocal pseudoptential operator for the case ofextended/condensed matter systems. The trace in eq. 78 is to be interpreted as being over trace-class operators on L ( C ) . Since ˜Γ is traceclass and V nl D is a finite rank operator (and hence bounded), their product is trace class as well, thusmaking eq. 78 well defined. By making use of the fact that the projection functions χ k ; p are related to atomic orbitals (and areusually expressed as the product of a spherical harmonic with a compactly supported radially symmetricfunction), we may move the action of Υ m h in eq. 81 from x , to the atomic positions x k and rewrite theabove in a more computationally convenient form as [33, 60]: ˆ χ k,p ( x ; η ; x k ) = (cid:88) m ∈ Z χ k ; p ( x ; Υ m h ◦ x k ) e − i πmη . (82)This modified form can be immediately verified for example when projection functions from s-orbitalsare involved since in that case χ k ; p ( x ; x k ) = χ k ; p ( (cid:107) x − x k (cid:107) R ) . K (Λ , Ψ , P G , D , G ) = 2 Tr. [ U V nl D U − Γ ]= 2
Tr. (cid:34)(cid:18) M G (cid:88) k =1 (cid:88) p ∈N k γ k,p ˆ χ k,p ( · ; η ; x k ) ˆ χ k,p ( · ; η ; x k ) (cid:19)(cid:18) (cid:90) ⊕ I ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1) ψ j ( x ; η ) ψ j ( y ; η ) dη (cid:19)(cid:35) . (83)Now using the definition of the trace and that the helical Bloch states ψ j ( x ; η ) are a basisof L ( D ) , this reduces to: K (Λ , Ψ , P G , D , G )= 2 M G (cid:88) k =1 (cid:88) p ∈N k γ k,p (cid:90) I (cid:18) ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1)(cid:10) ˆ χ k,p ( · ; η ; x k ) , ψ j ( · ; η ) (cid:11) L ( D ) (84) × (cid:10) ˆ χ k,p ( · ; η ; x k ) , ψ j ( · ; η ) (cid:11) L ( D ) (cid:19) dη = 2 M G (cid:88) k =1 (cid:88) p ∈N k ∞ (cid:88) j =1 γ k,p (cid:90) I (cid:18) f T e (cid:0) λ j ( η ) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:10) ˆ χ k,p ( · ; η ; x k ) , ψ j ( · ; η ) (cid:11) L ( D ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:19) dη = 2 M G (cid:88) k =1 (cid:88) p ∈N k ∞ (cid:88) j =1 γ k,p (cid:90) I (cid:18) g j ( η ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:10) ˆ χ k,p ( · ; η ; x k ) , ψ j ( · ; η ) (cid:11) L ( D ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:19) dη (85)Electrostatic Energy Term: We now discuss the contribution of the electrostatic in-teraction energy to the free energy per unit fundamental domain. This is the fourth termon the right-hand side of eq. 67. Often, it is computationally advantageous to expressthis term using a so-called local formulation [60, 113–115], as we now do . For a finitesystem with atomic nuclei located at the points { p k ∈ R } N at k =1 and electron density ρ ,this term takes the form of the following optimization problem in the total electrostaticpotential Φ : E finiteel = max Φ (cid:26) − π (cid:90) R |∇ Φ | d x + (cid:90) R ( ρ + b finite )Φ d x (cid:27) + E finitesc ( p , p , . . . , p k ) . (86)Here b finite represents the total nuclear pseudocharge for the finite set of nuclei, and canbe expressed in terms of the individual nuclear pseudocharges (cid:8) b k ( x ; p k ) (cid:9) N at k =1 as: b finite ( x ) = N at (cid:88) k =1 b k ( x ; p k ) , (87)and the term E finitesc ( p , p , . . . , p k ) corrects for self-interactions and overlaps of the nu-clear pseudocharges [114]. Note that by design, the individual nuclear pseudocharges are The term local formulation is associated with the fact that the electrostatic field Φ can be solvedthrough a Poisson equation, which avoids evaluation of the non-local integrals in eq. 88. Φ that solves the maximization problem in 86 is theNewtonian potential associated with the net charge in the system: Φ( x ) = (cid:90) R ρ ( y ) + b finite ( y ) (cid:107) x − y (cid:107) R d y (88)To extend the above formulation to a helical structure, we first write the total nuclearpseudocharge at any point x ∈ C in terms of the pseudocharges of the atoms in thefundamental domain as: b ( x , P G , G ) = (cid:88) m ∈ Z M G (cid:88) k =1 b k ( x ; Υ m h ◦ x k ) , (89)and observe that this quantity is group invariant owing to the aforementioned propertiesof the individual nuclear pseudocharges [39]. Since the electron density is group invariantas well, it follows that the net electrostatic potential Φ expressed as : Φ( x ) = (cid:90) C ρ ( y ) + b ( y , P G , G ) (cid:107) x − y (cid:107) R d y , (90)is also group invariant [8]. Therefore, we may use Φ , b and ρ as defined over the funda-mental domain, to define the analog of eq. 86 for the helical structure as: E el ( ρ, P G , D , G ) = max Φ (cid:26) − π (cid:90) D |∇ Φ | d x + (cid:90) D (cid:18) ρ ( x ) + b ( x , , P G , G ) (cid:19) Φ( x ) d x (cid:27) + E sc ( P G , G , D ) . (91)For the sake of brevity, we omit the details of the form of the corrections due to selfinteractions and overlaps of the nuclear pseudocharges, as reduced to the fundamentaldomain (i.e., the term E sc ( P G , G , D ) above) and instead point to [39, 114] for relevantinformation.Electronic Entropy Term: Finally, the last term on the right-hand side of 67 representsthe electronic entropy contribution to the free energy. Following [116], this term can beexpressed for a finite system with density matrix D as: S finite = − k B Tr. [ D log( D ) + ( I − D ) log( I − D )] , (92)with I denoting the identity operator. To obtain the analogous expression for the case ofthe helical structure, we work with the trace per unit fundamental domain instead. Thisgives us: S = − k B Tr. (cid:20) ˜Γ log(˜Γ) + (
I − ˜Γ) log(
I − ˜Γ) (cid:21) = − k B (cid:90) I Tr. [˜Γ η log(˜Γ η ) + ( I − ˜Γ η ) log( I − ˜Γ η )] dη . (93) Using a Fourier series expansion, it can be shown (see e.g. [43] for similar arguments) that eq. 90 iswell defined whenever the system is charge neutral i.e., when (cid:90) D ρ ( x ) + b ( x , P G , G ) d x = 0 .
26y means of spectral mapping [91], the use of eq. 66, and by noting again that the helicalBloch states ψ j ( x ; η ) are a basis of L ( D ) , we readily obtain: S (Λ) = − k B (cid:90) I (cid:20) ∞ (cid:88) j =1 f T e (cid:0) λ j ( η ) (cid:1) log (cid:0) f T e (cid:0) λ j ( η ) (cid:1)(cid:1) + (cid:0) − f T e (cid:0) λ j ( η ) (cid:1)(cid:1) log (cid:0) − f T e (cid:0) λ j ( η ) (cid:1)(cid:1)(cid:21) dη = − k B (cid:90) I (cid:20) ∞ (cid:88) j =1 g j ( η ) log (cid:0) g j ( η ) (cid:1) + (cid:0) − g j ( η ) (cid:1) log (cid:0) − g j ( η ) (cid:1)(cid:21) dη (94)With the above terms explicitly defined, we now turn to the variational problem fordeducing the governing equations.Variational Problem and Kohn-Sham Equations: The variational problem for Kohn-Sham ground state of a given helical structure (i.e., the atomic coordinates in the fun-damental domain, the fundamental domain and the helical group are held fixed) con-sists of minimizing the electronic free energy F (Λ , Ψ , P G , D , G ) with respect to thehelical Bloch states and the helical bands, subject to the constraint in eq. 64. In theliterature, this minimization is often stated in terms of the helical Bloch states andthe electronic occupation numbers instead [33, 60]. Along those lines, we may define (cid:101) F ( G , Ψ , P G , D , G ) = F (Λ , Ψ , P G , D , G ) to write the variational problem as: (cid:101) F ( P G , D , G ) = inf. Ψ , G (cid:101) F ( G , Ψ , P G , D , G ) (95)subject to: (cid:90) I ∞ (cid:88) j =1 g j ( η ) = N e , (96)and the requirement that the states in Ψ be helical Bloch states. This requires that forany ψ j ( · ; η ) ∈ Ψ and m ∈ Z , we have ψ j (Υ m h ◦ x ; η ) = e − i πmη ψ j ( x ; η ) , as well as theorthonormality condition between two helical Bloch states ψ i ( · ; η ) , ψ j ( · ; η ) ∈ Ψ : (cid:104) ψ i ( · ; η ) , ψ j ( · ; η ) (cid:105) L ( D ) = δ i,j . (97)We take variations of the above constrained minimization problem, and obtain the Euler-Lagrange equations as the following helical symmetry adapted Kohn-Sham equations overthe fundamental domain: H KS [ G , Ψ , P G , D , G ] ψ j ( · ; η ) = λ j ( η ) ψ j ( · ; η ) for j ∈ N , η ∈ I . (98)Here, the helical symmetry adapted Kohn-Sham Hamiltonian operator (with its depen-dence on the helical Bloch sates, the occupation numbers, etc., made explicit) is: H KS [ G , Ψ , P G , D , G ] ≡ −
12 ∆ + V xc + Φ + (cid:101) V nl D , (99) Up to a notational change, this operator is essentially the same as H KS [Λ , Ψ , P G , D , G ] , when the depen-dence on helical bands (instead of the occupation numbers) is highlighted.
27n which V xc = δE xc ( ρ, D ) δρ is the exchange correlation potential, Φ is the net electro-static potential and satisfies the following symmetry adapted Poisson problem over thefundamental domain: − ∆Φ = 4 π (cid:0) ρ + b ( · , P G , G ) (cid:1) , Φ(Υ h ◦ x ) = Φ( x ) , (100)and the operator (cid:101) V nl D = U V nl D U − is as defined in eqs. 80 and 81.Harris-Foulkes Functional: The above set of expressions represent a set of coupled non-linear partial differential equations in the fields ψ j ( · ; η ) and the scalars g j ( η ) . Once theyhave been solved self-consistently, the ground state electronic free energy (cid:101) F ( P G , D , G ) per unit fundamental domain can be computed through eq. 67. In practical calcula-tions, since self-consistency is never achieved perfectly, a better estimate of the groundstate electronic free energy may be found using the so-called Harris-Foulkes functional[117, 118]. This can be written in helical symmetry-adapted form, using quantities ex-pressed over the fundamental domain as: F HF (Λ , Ψ , P G , D , G ) = E band (Λ) + E xc ( ρ, D ) − (cid:90) D V xc ( ρ ( x )) ρ ( x ) d x + 12 (cid:90) D (cid:18) b ( x , , P G , G ) − ρ ( x ) (cid:19) Φ( x ) d x + E sc ( P G , G , D ) − T e S (Λ) . (101)All the quantities on the right-hand side of the above equation are easily interpreted basedon earlier discussion, except the first one, i.e., E band (Λ) , which represents the electronicband energy per unit fundamental domain. For a finite system with a single electronHamiltonian H and single particle density matrix D , this quantity is expressed as [103]: E finiteband = 2 Tr. [ H D ] . (102)Analogously, for the helical structure, we use the trace per unit fundamental domain towrite: E band = 2 Tr. [ H ˜Γ ] = 2 (cid:90) I Tr. [ H η ˜Γ η ] dη . (103)Using eqs. 58 and 66 and using the completeness of the helical Bloch waves, we see thatthis is be expressible as: E band (Λ) = 2 (cid:90) I ∞ (cid:88) j =1 λ j ( η ) f T e (cid:0) λ j ( η ) (cid:1) dη = 2 (cid:90) I ∞ (cid:88) j =1 λ j ( η ) g j ( η ) dη (104)Atomic Forces: The Hellman-Feynman forces on the atoms in the fundamental domainare (in Cartesian coordinates): f k = − ∂ (cid:101) F ( P G , D , G ) ∂ x k (105)28y directly differentiating the various terms involved (eqs. 83, 91) we arrive at the fol-lowing expression for f k using quantities specified over the fundamental domain : f k = (cid:88) m ∈ Z ( R πmα ) − (cid:90) D ∇ b k ( x ; Υ m h ◦ x k )Φ( x ) d x − ∂E sc ( P G , G , D ) ∂ x k − ∞ (cid:88) j =1 (cid:32) (cid:90) I g j ( η ) (cid:88) p ∈N k γ k ; p Re. (cid:40)(cid:20) (cid:90) D ˆ χ k,p ( x ; η ; x k ) ψ j ( x ; η ) d x (cid:21) × (cid:20) (cid:90) D ψ j ( x ; η ) ∂ ˆ χ k,p ( x ; η ; x k ) ∂ x k d x (cid:21)(cid:41)(cid:33) dη (107)Here, Re. denotes the real part of the quantity in braces.This completes a discussion of the derivation of the various physically relevant terms,as well as the form of the equations of Kohn-Sham theory, as applied to a helical structureassociated with a helical group generated by a single element. Comments on modificationsof the above equations while dealing with a structure associated with a helical groupgenerated two elements, and a presentation of the final expressions/equations for thatcase in appear in Appendix C.
3. Numerical Implementation
The Kohn-Sham equations for a helical structure (i.e., eq. 98 or eq. C.15) are a set ofnon-linear eigenvalue problems indexed by η (as well as ν in case of the group G ) thatare coupled to each other through the electron density ρ . The standard procedure forsolving the equations of Kohn-Sham theory is through self-consistent field (SCF) itera-tions [24]. This amounts to starting from a reasonable guess of the electron density in thefundamental domain (e.g. superpositions of individual atomic densities, as is used in oursimulations) and an appropriate set of trial orthonormal wavefunctions (randomly chosenin our simulations), and then evaluating the eigenstates of the Kohn-Sham operator withthese guesses. Thus, a set of linear eigenvalue problems (i.e., those associated with thelinearized Kohn-Sham operator evaluated at the given electron density) indexed by η (aswell as ν in case of the group G ) have to be solved. From this, the (trial) Fermi-level ofthe system and the (trial) occupation numbers maybe computed. The eigenfunctions andthe occupation numbers may be then combined (in accordance with eq. 61 or eq. C.3) toyield the trial electron density for the next step of the iterations. The above procedurehas to be repeated till the difference in the electron density (or the effective potential)between successive iterations reaches below the desired convergence threshold. We will Motivated by [60, 119] we may use integration by parts to modify the last term on the right-hand sideof eq. 107, so that the derivatives of the projectors with respect to atomic coordinates can be eliminatedin favor of Cartesian gradients of the wavefunctions instead. This tends to improve the accuracy of thecomputed forces in practical calculations – the orbitals are more smoothly varying than the projectorsand therefore they tend to behave better upon taking derivatives. With this change as well as makinguse of the discussion in Footnote 34, we have: (cid:90) D ψ j ( x ; η ) ∂ ˆ χ k,p ( x ; η ; x k ) ∂ x k d x = (cid:88) m ∈ Z ( R πmα ) − (cid:90) D ∇ ψ j ( x ; η ) ˆ χ k,p ( x ; η ; Υ m h ◦ x k ) e i πmη d x . (106) Many quantities described in Section 2.2.2 and Appendix C involve integrals over η ∈ I = [ − , ) (as well as normalized summations over ν ∈ { , , , . . . , N − } for thegroup G ). To evaluate such integrals numerically, we employ quadrature based on theMonkhorst-Pack scheme [120]. Specifically, we sample the interval I using a grid of N η points, and write: (cid:90) I f ( η ) dη ≈ N η (cid:88) b =1 w b f ( η b ) . (108)Here, w b and η b denote the integration weights and integration nodes respectively. Sum-mations over ν are left unchanged. The total number of points used for discretizing thereciprocal space (i.e., set B = I × { , , , . . . , N − } ), therefore, is N K = N η × N (with N = 1 for the group G ). Based on considerations of time-reversal symmetry (whichapply as long as e.g. magnetic fields are absent) [60, 121], it follows that for η ∈ I and ν = 1 , , N − : λ j ( η, ν ) = λ j ( − η, N − ν ) , ψ j ( x ; η, ν ) = ψ j ( x ; − η, N − ν ) , (109)while for ν = 0 : λ j ( η,
0) = λ j ( − η, , ψ j ( x ; η,
0) = ψ j ( x ; − η, . (110)Effectively, the above considerations reduce the number of quadrature points over recip-rocal space by half (i.e., N K ≈ N η × N ).With the above discretization choices, the self-consistent field iterations for the Kohn-Sham problem amount to solving a series of N K linear eigenvalue problems on everyiteration step. Based on the mathematical treatment presented earlier as well as in[39], it follows that eigenvalue problems associated with distinct values of η (i.e. η b indiscretized form) and/or ν are disjoint from each other. This implies that these (linear)eigenvalue problems can be solved independently of each other, in an embarrassinglyparallel manner, regardless of how the discretization in real space is carried out. Wemake use of this feature of the equations to assign these distinct eigenvalue problems todifferent computational cores. This serves as a natural parallelization scheme and helpsin drastically reducing the wall time associated with the most computationally intensivepart of the SCF iterations. Primarily, there are two distinct sources of infinite sums in the equations presented inSection 2.2.2 and Appendix C. The first arises due to summing over an infinite number ofhelical bands (e.g., eqs. 61 and C.3). To truncate such sums we assume that the electronicoccupation numbers reduce to zero beyond the lowest N s bands and therefore, only N s eigenstates for each value of η b (and also each ν for G ) need to be computed duringthe self-consistent field iterations. In effect, this is also an enforcement of the Aufbau30rinciple for the system [25, 43]. Depending on the size of the discretized reciprocal space(i.e., the value of the number N K ) we have found that including just a few extra bandsbeyond the minimum number required for holding the N e electrons per unit fundamentaldomain, suffices. The second source of infinite sums arises from considering terms associated with grouporbits (e.g., eqs. 89 and C.12), since helical groups by definition are infinite. However,these sums are also always associated with functions that are supported in a small ballaround an atom of the structure (e.g. the nuclear pseudocharge in eq. 89 and the nonlocalpseudopotential projection function in eq. 81). Therefore, the influence of such sums onpoints in the fundamental domain is only dependent on terms of the summation thatresult in a nonzero overlap between the function support and the fundamental domain.This allows such infinite sums to be truncated as well.
To carry out a discretization of the governing equations in a manner that is naturallyadapted to the underlying helical symmetries of the structures being studied, it is usefulto employ helical coordinates, as introduced in [8]. In order to have this coordinatesystem be commensurate with the helical groups G or G , the coordinate transformationformulae are as follows. For x ∈ C , if the Cartesian coordinates are ( x , x , x ) , then thecorresponding helical coordinates ( r, θ , θ ) are: r = (cid:113) x + x , θ = x τ , θ = 12 π arctan (cid:18) x x (cid:19) − α x τ . (111)The coordinates ( r, θ , θ ) are a natural generalization of cylindrical coordinates in thesense that they effectively reduce to cylindrical coordinates when the twist angle parame-ter α of the system is set to zero. We may verify that these relations are onto and globallyinvertible on C\{ t e : t ∈ R } . Furthermore, the inverse relations: ( r, θ , θ ) (cid:55)→ ( x , x , x ) = (cid:0) r cos(2 π ( αθ + θ )) , r sin(2 π ( αθ + θ )) , θ τ (cid:1) (112)map the open cuboid (0 , R ) × (0 , × (0 , to the interior of the fundamental domain of G , (i.e., to the set D ) and the open cuboid (0 , R ) × (0 , × (0 , / N ) to the interior ofthe fundamental domain of G (i.e., to the set (cid:101) D ).The action of a helical group can be easily computed in these coordinates as follows.Let x = ( x , x , x ) ∈ C have helical coordinates ( r, θ , θ ) . The action of the isometry Υ h that generates G (and also G ) is to map x to the point x (cid:48) = Υ h ◦ x = R πα x + τ e . Denoting the Cartesian and helical coordinates of this new point as ( x (cid:48) , x (cid:48) , x (cid:48) ) and ( r (cid:48) , θ (cid:48) , θ (cid:48) ) , respectively, we see that x (cid:48) = x cos(2 πα ) − x sin(2 πα ) , x (cid:48) = x cos(2 πα ) + x sin(2 πα ) and x (cid:48) = x + τ . Now, using eqs. 112 and 111, we get r (cid:48) = r, θ (cid:48) = θ + 1 , θ (cid:48) = θ . By a similar calculation, we see that the action of the second generator Υ c of the This is a well used approximation strategy in the literature (see e.g. [32, 33, 36]). For finite systemsat electronic temperatures that are less than a few thousand Kelvin, it suffices to choose the numberof states to be equal to a multiple of half the number of electrons, with the multiplication factor beingbetween . and . [122]. For an extended system like a helical structure, often a just a few extrabands beyond half the number of electrons is sufficient since this actually amounts to these few extrastates being available for every value of η b or ν . a]As a result of this, tens or even hundreds of extrastates (with occupation numbers approaching zero) get effectively included in the calculations. G (i.e., the pure rotation by π N about axis e ), is to map the point with helicalcoordinates ( r, θ , θ ) to the point ( r, θ , θ + N ) . These calculations imply in particularthat if a function is invariant under the group G , then it is periodic in θ , with period ,when expressed in helical coordinates. Similarly, invariance of a function under the group G , implies periodicity in θ (with period ), as well as periodicity in θ , (with period / N )when expressed in helical coordinates. These observations make it easy to enforce helicalBloch boundary conditions on the wavefunctions, as well as the group invariance of theelectrostatic potential, in simulations. Derivation of the Cartesian gradient operator, theLaplacian operator and the volume integral in helical coordinates appears in AppendixD. We employ a finite difference strategy for discretizing the governing equations in realspace. We choose an annular cylindrical region Ω with axis along e as the simulationdomain. This allows us to avoid the singularity associated with the helical coordinatesystem along the axis e and does not present any issues as long as the atoms of thestructure are located well within the annular region [39, 60]. This latter condition is wellsatisfied by the nanotube structures simulated in this work. The set Ω can also serveadequately as a suitable fundamental domain for either group G or G in simulations(i.e., it can replace D or (cid:101) D in the formulae presented in Section 2.2.2 and Appendix C).In cylindrical coordinates ( r, ϑ, z ) we have: Ω = (cid:8) ( r, ϑ, z ) ∈ R : R in ≤ r ≤ R out , ≤ ϑ ≤ Θ , ≤ z ≤ τ (cid:9) . (113)The boundary of Ω can be expressed as: ∂ Ω = ∂R in (cid:91) ∂R out (cid:91) ∂ϑ (cid:91) ∂ϑ Θ (cid:91) ∂ Z (cid:91) ∂ Z τ , (114)with ∂R in and ∂R out denoting the surfaces r = R in and r = R out respectively, ∂ϑ and ∂ϑ Θ denoting the surfaces ϑ = 0 and ϑ = Θ respectively, and finally, ∂ Z and ∂ Z τ denoting the surfaces z = 0 and z = τ respectively. Figure 3 illustrates the simulationcell as well as the boundaries of this domain.We set up a finite difference grid over Ω using helical coordinates, with spacing h r , h θ and h θ along the r , θ and θ directions respectively. Since the simulation domaincan be represented in helical coordinates with r ranging from R in to R out , θ rangingfrom to , and θ ranging from to N , it follows that R out − R in = N r h r , τ = N θ h θ ,and N = N θ h θ , for natural numbers N r , N θ and N θ . We index the finite differencenodes using a triplet of natural numbers ( i , j , k ) , for i = 1 , , . . . , N r , j = 1 , , . . . , N θ and k = 1 , , . . . , N θ . We denote the value of a function f at the grid point ( i , j , k ) as f ( i , j , k ) .We will denote h = Max. (cid:18) h r , τ h θ , (cid:0) R in + R out (cid:1) h θ (cid:19) as the mesh spacing.The formulae presented in Section 2.2.2 and Appendix C require Cartesian gradients,the Laplacian operator and integrals over Ω to be evaluated using the finite differencescheme. The formulae for these quantities, as expressed in helical coordinates appearsin Appendix D. With these in hand, we approximate the first order partial derivatives32 XZ τ R out R in Θ ∂ Z τ ∂ Z ∂ ϑ Θ ∂ R in ∂ R out ∂ ϑ Figure 3: Illustration of the simulation domain Ω (domain boundary lines in blue) containing a fewatoms, as well as the various bounding surfaces of the domain. For a structure associated with a two-generator helical group, the parameter Θ = 2 π/ N relates to the cyclic group order, while τ is related tothe pitch of the screw transformation. using central differences as: ∂f∂r (cid:12)(cid:12)(cid:12)(cid:12) ( i , j , k ) ≈ n o (cid:88) p =1 (cid:18) w first p,r (cid:0) f ( i + p, j , k ) − f ( i − p, j , k ) (cid:1)(cid:19) ,∂f∂θ (cid:12)(cid:12)(cid:12)(cid:12) ( i , j , k ) ≈ n o (cid:88) p =1 (cid:18) w first p,θ (cid:0) f ( i , j + p, k ) − f ( i , j − p, k ) (cid:1)(cid:19) ,∂f∂θ (cid:12)(cid:12)(cid:12)(cid:12) ( i , j , k ) ≈ n o (cid:88) p =1 (cid:18) w first p,θ (cid:0) f ( i , j , k + p ) − f ( i , j , k − p ) (cid:1)(cid:19) . (115)We approximate the pure (i.e., non-mixed) second order partial derivatives as: ∂ f∂r (cid:12)(cid:12)(cid:12)(cid:12) ( i , j , k ) ≈ n o (cid:88) p =0 (cid:18) w second p,r (cid:0) f ( i + p, j , k ) + f ( i − p, j , k ) (cid:1)(cid:19) ,∂ f∂θ (cid:12)(cid:12)(cid:12)(cid:12) ( i , j , k ) ≈ n o (cid:88) p =0 (cid:18) w second p,θ (cid:0) f ( i , j + p, k ) + f ( i , j − p, k ) (cid:1)(cid:19) ,∂ f∂θ (cid:12)(cid:12)(cid:12)(cid:12) ( i , j , k ) ≈ n o (cid:88) p =0 (cid:18) w second p,θ (cid:0) f ( i , j , k + p ) + f ( i , j , k − p ) (cid:1)(cid:19) , (116)while the mixed second order partial derivative ∂ f∂θ ∂θ is obtained by applying the above33rst order derivative formula first in θ and then in θ , i.e., ∂ f∂θ ∂θ (cid:12)(cid:12)(cid:12)(cid:12) ( i , j , k ) ≈ n o (cid:88) p =1 w first p,θ (cid:20)(cid:18) ∂f∂θ (cid:19) ( i , j , k + p ) − (cid:18) ∂f∂θ (cid:19) ( i , j , k − p ) (cid:21) ≈ n o (cid:88) p =1 w first p,θ (cid:20)(cid:26) n o (cid:88) p (cid:48) =1 w first p (cid:48) ,θ (cid:0) f ( i , j + p (cid:48) , k + p ) − f ( i , j − p (cid:48) , k + p ) (cid:1)(cid:27) − (cid:26) n o (cid:88) p (cid:48) =1 w first p (cid:48) ,θ (cid:0) f ( i , j + p (cid:48) , k − p ) − f ( i , j − p (cid:48) , k − p ) (cid:1)(cid:27)(cid:21) (117)In the above formulae, n o denotes half the finite difference order and is set to for allour simulations (i.e., th order finite differences are used). This choice has also beenemployed elsewhere [32, 33, 39, 60] and is found to be adequate for attaining chemicalaccuracy. Letting s denote r , θ or θ , the weights that appear in the above formulae canbe expressed as [123]: w second ,s = − h s n o (cid:88) q =1 q ,w second p,s = 2( − p +1 h s p ( n o !) ( n o − p )!( n o + p )! for p = 1 , , . . . , n o ,w first p,s = ( − p +1 h s p ( n o !) ( n o − p )!( n o + p )! for p = 1 , , . . . , n o . (118)We employ the following quadrature rule for approximating integrals over Ω : (cid:90) Ω f ( x ) d x ≈ h r h θ h θ N r (cid:88) i =1 N θ (cid:88) j =1 N θ (cid:88) k =1 πτ r i f ( i , j , k ) , (119)with r i denoting the radial coordinate of the finite difference node ( i , j , k ) . We need to specify the boundary conditions on the various fields that are being solvedfor in the governing equations. These are the helical Bloch states ψ j ( x ; η, ν ) (for η ∈ I and ν = { , , , . . . , N − } ) which satisfy the Kohn-Sham equations (eq. 98 or eq. C.15), andthe total electrostatic potential Φ which satisfies Poisson’s equation (eq. 100 or eq. C.20).In helical coordinates, we may interpret the helical Bloch boundary conditions (eq. C.1)as the following conditions: ψ j ( r, θ = 1 , θ ; η, ν ) = e − i πη ψ j ( r, θ = 0 , θ ; η, ν ) , (120)which applies to the surfaces ∂Z (cid:83) ∂Z τ ; as well as the condition: ψ j ( r, θ , θ = 1 N ; η, ν ) = e − i π ν N ψ j ( r, θ , θ = 0; η, ν ) , (121)which applies to the surfaces ∂ϑ (cid:83) ∂ϑ Θ . We assume that the atoms within Ω are suffi-ciently far away from the boundary surfaces ∂R in and ∂R out , so that the electron density34ecays to zero at these surfaces and zero Dirichlet boundary conditions on the wavefunc-tions (Section 2.2.1) can be applied. Thus, the conditions to be applied on the surfaces ∂R in (cid:83) ∂R out are: ψ j ( r = R in , θ , θ ; η, ν ) = ψ j ( r = R out , θ , θ ; η, ν ) = 0 . (122)As already discussed, the electrostatic potential Φ inherits the symmetry of the he-lical structure (i.e., it is invariant under G or G ). Thus, it follows that the boundaryconditions on this quantity on the surfaces ∂Z (cid:83) ∂Z τ and ∂ϑ (cid:83) ∂ϑ Θ are respectively: Φ( r, θ = 1 , θ ) = Φ( r, θ = 0 , θ ) , Φ( r, θ , θ = 1 N ) = Φ( r, θ , θ = 0) . (123)To apply the right boundary conditions on Φ on the surfaces ∂R in (cid:83) ∂R out , we mayevaluate eq. 88 or eq. C.13 directly using Ewald summation [71] or multipole expansion[88] techniques. On every SCF iteration step, the lowest N s eigenstates of a set of N K Kohn-Shamoperators (indexed by η b for the case of G , and η b , ν for the case of G ) have to becomputed. We employ iterative diagonalization based on Chebyshev polynomial filteredsubspace iterations (CheFSI) [124–126] for this purpose. Due to the fact that the helicalcoordinate system is curvilinear, the finite difference discretization of the Laplacian op-erator results in a discretized operator that is non-Hermitian (even though the Laplacianas an operator on an infinite dimensional Hilbert space is Hermitian). This results in thediscretized Hamiltonian operator also being non-Hermitian, which can lead to non-realeigenvalues of the matrix. However, as the discretization is made finer (i.e., h r , h θ , h θ become smaller) the discretized Laplacian and Hamiltonian matrices approach Hermi-tian matrices and so, the eigenvalues resulting from these discretized operators tend tohave vanishingly small imaginary parts [39, 127]. In practice, for the mesh spacings thathave been considered in this work, the imaginary parts are small enough that they canbe safely ignored without affecting the quality of the simulations (also see [39] and [60]where a similar situation was encountered).We use Chebyshev polynomial filter orders in the range to for our simulations.This is somewhat higher than what is commonly employed in finite difference DFT cal-culations in affine coordinate systems with comparable mesh spacing [32, 33, 125]. Weadopt it here to mitigate the effect of the larger spectral width of the discretized Hamil-tonian that arises due to crowding of grid points as one approaches the origin in helicalcoordinates. We have also observed that the time for computing matrix-vector productsusing the discretized Hamiltonian in helical coordinates is larger compared to the timerequired for the case of a finite difference Hamitonian (of the same size) arising froma cylindrical system (obtained by setting the twist angle parameter α to in the heli-cal case). This is certainly due to the presence of cross derivatives in the helical case As remarked in [32], it is often adequate to set the net potential to zero on the bounding surfaces,instead of applying the above procedures, provided the boundaries are sufficiently distant. The netpotential arises from a net neutral charge distribution, and therefore it usually decays to zero faster.This is the procedure adopted here for the simulations presented in Section 4. Φ . To accelerate convergence, we use anincomplete LU factorization based preconditioner [129]. We have implemented the above computational strategies using the MATLAB [130]software package into a code called Helical DFT. The code parallelizes computation overthe different η and ν values using MATLAB’s Parallel Computing Toolbox (the parfor function). For maximum efficiency of the MATLAB implementation, code vectorizationhas been used as much as possible. However, for computing certain quantities (suchas the atomic forces and the net nuclear pseudocharge), multiple levels of nested loopswere found unavoidable. These routines were converted into machine code by use ofthe MATLAB Coder framework, which helped alleviate performance issues. In orderto reduce the memory footprint associated with the storage of the different Hamiltonianmatrices arising from the N K different values of η b (and also ν for G ), we avoid computingmatrix vector products through MATLAB’s internal sparse matrix framework since thatrequires these matrices to be available explicitly. Instead, we store only the Laplacianpart of the Hamiltonian matrix in compressed sparse row (CSR) format and apply thehelical Bloch boundary conditions associated with the different values of η b and/or ν on the fly, while computing matrix vector products. The actual task of computing thesematrix vector products is carried out using a C language routine which has been compiledand interfaced with our MATLAB code.Some other relevant details related to the implementation are as follows. We usethe periodic variant of Pulay’s scheme [131, 132] in the total potential to accelerate theconvergence of the SCF iterations. The Fermi energy is calculated using a nonlinearequation root finder (MATLAB’s fzero function). When required, structural relaxation isachieved by an implementation of the Fast Intertial Relaxation Engine (FIRE) algorithm[133].
4. Simulation Results and Discussion
We now turn to a discussion of numerical simulations and results. All simulations wererun using a single node of the Mesabi cluster at the Minnesota Supercomputing Institute,or a single node of the Hoffman2 cluster at UCLA’s Institute for Digital Research andEducation. Each compute node of Mesabi has Intel Haswell E5-2680v3 processorsoperating at . GHz, and GB to TB of RAM. Each compute node of the Hoffman2cluster has two -core Intel Xeon Gold 6140 processors (with . MB cache, andrunning at . GHz), and
GB of RAM. Earlier versions of this MATLAB code were developed in collaboration with Phanish Suryanaryana,Georgia Institute of Technology. Details on a more efficient C/C++ implementation by Suryanarayanaand collaborators appears in their forthcoming work. T e = 315 . Kelvin was used for Fermi-Diracsmearing to help accelerate SCF convergence.The large majority of the simulations here have focused on the study of single wallnanotubes. Starting from the sheet of an elemental two-dimensional material, nanotubesof any chirality can be formed using the so-called “roll-up” construction [135]. Thisprocedure also allows us to see [7, 48] that such nanotubes can be adequately representedusing helical groups generated by two elements with just atoms in the fundamentaldomain [60]. In this representation, the twist angle parameter α becomes related to thechirality of the tubes ( α = 0 for achiral tubes), while the parameter N , associated withthe cyclic group order, is related to the tube radius [48]. In our nanotube simulations,we have ensured that the atoms within the fundamental domain are always located to Bohrs away from the boundary surfaces ∂R in and ∂R out so as to allow sufficient decayof the electron density and the wavefunctions in the radial direction. Single-layer black phosphorus, or phosporene, is a two-dimensional nanomaterial thathas been the object of intense investigation in recent years due its association with a num-ber of unusual and fascinating material properties [136–142]. Nanotubes of this material,as formed by the roll-up construction have also received attention in the literature [143–156], due to their interesting optical and electronic properties, and the coupling of theseproperties to mechanical strains. This motivates our choice in selecting this material forthe simulations presented in this work.As a starting point, we obtained the ground state structure of a single layer of phos-phorene (Figure 4) using the same pseudopotential, exchange correlation functional andelectronic temperature as the Helical DFT simulations subsequently described. We usedthe plane-wave DFT code ABINIT [29, 157] to perform the geometry relaxation calcula-tion. The periodic unit cell for this simulation contained atoms. An energy cutoff of Ha, along with × × k-points and a cell vacuum of Bohr in the Z-directionwas used. At the end of this very refined calculation the atomic forces were all less than − Ha/Bohr, while the cell stresses were of the order of − Ha/Bohr or lower. Someof the structural parameters obtained from the calculation are shown in Table 1. Thereappears to be generally good agreement with the literature thus giving us confidence inthe reliability of the subsequent simulations with regard to materials physics.To represent the zigzag phosphorene nanotubes in the Helical DFT simulations de-scribed next, we roll up the phosphorene sheet along the X-axis and place the atoms fromthe aforementioned periodic unit cell (the shaded region in Figure 4) into the Helical DFTsimulation cell ( Ω ). In the absence of relaxation effects, the pitch τ associated with thetwo-generator helical group, is equal to the lattice constant along the Y axis. Further-more, the angle Θ associated with the cyclic group order N is related to the radius of thenanotube via the relation Θ = 2 π N = a R avg. . Here a denotes the lattice vector along the The minor differences are possibly due to our use of LDA exchange correlation, which tends topredict overbinding and shortened bond lengths [158]. attice constant a Lattice constant a Bond length δ Bond length δ Bond angle γ along X axis (Å) along Y axis (Å) (Å) (Å) (degrees) .
26 4 .
36 2 .
20 2 .
193 102 . ( . [139], . [159]) ( . [139], . [159]) ( . [159]) ( . [159]) ( . [159]) Table 1: Structural parameters of phosphorene computed using a periodic DFT calculation. Quantitiesin parentheses refer to values in the literature along with references. Explanation of the parameters isavailable from Figure 4 . XY Simulation cell a a τ = a (a) Top view of a phosphorene sheet. Atomsin the shaded region are placed in the HelicalDFT simulation cell (fundamental domain). δ δ γ YZX (b) Side view of a phosphorene sheet. Differ-ent structural parameters are labeled. Valuesavailable in Table 1.Figure 4: Phosphorene lattice as computed using a periodic DFT calculation. This provides the startingpoint for Helical DFT calculations involving phosphorene nanotubes. For zigzag tubes, the parameter τ in the two-generator helical group is set to be equal to the lattice constant a along the Y-axis. X-axis in the phosphorene sheet, and R avg. denotes the average radial coordinate of theatoms in the fundamental domain of the nanotube. For chiral nanotubes, we additionallyinclude a no-zero twist angle parameter α . Figure 5 shows examples of two phosphorenenanotubes studied in this work using Helical DFT. To study the convergence properties of our numerical implementation, we considera chiral phosphorene nanotube of radius . nanometers. The cyclic group order is N = 32 and the twist angle parameter is α = 0 . . For a given structure and a fixedsimulation domain, the two convergence parameters under study are the mesh spacing h = Max. (cid:18) h r , τ h θ , (cid:0) R in + R out (cid:1) h θ (cid:19) , and the number of points N η used for discretizingthe set I in reciprocal space. As a reference calculation, we computed the electronicstructure of the system with the finest mesh ( h = 0 . Bohr) and the highest value of N η ( = 23 ) that we could afford under computational resource constraints. First, forthe mesh convergence study, we fix N η = 23 and carry out a series of calculations with h ≈ . , . , . , . , . , . Bohr. Next, for studying convergence with respect to N η , we fix h = 0 . Bohr and carry out a series of calculations with N η = 1 , , , , , .We plot the errors in the energy per atom and the atomic forces in Figure 6 in each ofthe above cases in Figure 6.From the figures, it is clear that the code converges to the reference calculations38 YZ (a) Zigzag phosphorene nanotube. X YZ (b) Chiral phosphorene nanotube obtained by set-ting the angle of twist parameter α = 0 . .Figure 5: Representative phosphorene nanotubes (radius = . nanometers) studied in this work usingHelical DFT. Insets show zoomed-in views of the atomic arrangement to highlight differences betweenthe two cases. systematically. Using straight-line fits to the convergence data with respect to h , we findslopes of . and . for the energies and forces, respectively. These numbers are verynearly identical to the convergence rates observed in finite difference calculations involvingcylindrical coordinates [60] and are also comparable to finite difference calculations inaffine coordinate systems [32, 33]. From the data, we are also able to estimate thatthe parameters h = 0 . Bohr and N η = 11 are more than sufficient to reach chemicallyaccurate energies and forces (the exact errors for these choices with respect to the referencecalculation were × − Ha/atom and × − Ha/Bohr in the energies and forces,respectively). For computational efficiency therefore, we use this set of parameters forall relaxation calculations described subsequently, and switch to the reference calculationparameters (i.e., h = 0 . Bohr, N η = 23 ) only when more accurate energies / band gapsare required at the end of a relaxation procedure. Verifying the accuracy of the Helical DFT code with respect to standard plane-wavecodes can be challenging since the plane-wave codes may be required to include an enor-mous number of atoms in the periodic unit cell in order to replicate the exact system beingstudied by Helical DFT. In case of the above chiral nanotube system for example, over , atoms would be needed, making the planewave calculation unfeasible. Therefore,we carry out this accuracy check in two steps. First, we set up a Helical DFT calculationfor a zigzag nanotube (i.e., α = 0 ) of radius . nanometers. For this tube, the cyclicgroup order N = 16 . Then, we simulate this tube using the ABINIT code by employing The convergence thresholds (in terms of relative residuals) for the SCF iterations and the Poissonproblem had been set to − and × − for the above calculations, and we continue to use thesevalues for all subsequent simulations. atom unit cell (periodicity enforced along the Z axis, Dirichlet boundary conditionsalong X and Y axes by padding with a large amount of vacuum). We converged bothcodes to the extent allowed by computational resources, and observed that the energies(in Ha/atom) and the forces (in Ha/Bohr) from these two calculations agreed with eachother to × − or better. Next to study a case for which α (cid:54) = 0 , we set up an artificialsystem consisting of atoms along a single helix (similar to the configuration in Figure2). This system was generated using a single generator helical group ( N = 1 ) and used α = 0 . . The helical unit cell had atoms, while the periodic unit cell in ABINIT in-volved . Once again, upon convergence with respect to their respective discretizationparameters, the codes produced results that differed from each other by about × − Ha/atom or Ha/Bohr. This completes the accuracy tests. Through the above examples, we were also able to observe that for realistic helicalnanostructure simulations, the wall time for Helical DFT can be up to orders of magnitudesmaller compared to a well optimized plane-wave code like ABINIT, making it a powerfulfirst principles tool in the study of such systems. − − − − Mesh spacing h E rr o r i n e n e r g y a nd f o r ce s Error in energy (Ha/atom)Straight line fitError in forces (Ha/Bohr) (a) Convergence with respect to discretizationin real space. − − − − N η E rr o r i n e n e r g y a nd f o r ce s Error in energy (Ha/atom)Straight line fitError in forces (Ha/Bohr) (b) Convergence with respect to number ofpoints used to discretize I .Figure 6: Convergence properties of the Helical DFT code with respect mesh size h and the number ofpoints N η used to discretize the set I . The error in the forces is computed by considering the magnitudeof the maximum difference in all the force components of all the atoms. By using the two-generator helical group G , it is possible to describe torsional defor-mations in nanostructures [48, 52, 160–162]. We use this procedure here to illustrate It is likely that Helical DFT can be made to agree with ABINIT results to finer levels of accuracy.However, the quasi-one-dimensional nature of the systems being studied, and the slow convergence ofthe electrostatics requires the use of large amounts of vacuum padding in the ABINIT supercell, and thistends to cause serious convergence issues as the energy cutoff is increased. Due to the use of helical symmetries, the structures being modeled in the simulations are infinite. Inpractice, nanotube structures have finite extent, though some can have lengths of the order of macroscopicsizes [163]. The edge effects in these materials are expected to decay as one moves towards the interior ofthe material, both in the continuum elasticity sense [164, 165] and at the level of the electronic structure[166]. This offers a justification for the conceptual correctness of the simulations. ab initio .Specifically, we investigate the behavior of a zigzag phosphorene nanotube (nanotuberadius . nanometer, cyclic group order N = 64 .) under twisting deformations. Webegin with the untwisted structure ( α = 0 ) and relax the atoms in the simulation cell tillall force components on all atoms are below − Ha/Bohr. Starting from this relaxedstructure, we apply twisting deformations to the nanotube by prescribing non-zero valuesof the twist angle parameter α in the two-generator helical group G used for describingthe nanotubes. We varied α in steps of . and relaxed the resulting nanotube struc-ture in each case, till all components of the forces on all the atoms in the simulation celldropped below − Ha/Bohr once again. The largest twist we considered is about degrees per nanometer length of the tube, which corresponds to α = 0 . . Anticipatinga quadratic dependence of the nanotube twist energy for small values of alpha [48], weplot the results as shown in Figure 7. In that figure, we have also plotted the twistingenergies obtained when the atomic relaxation effects are not considered.We write the twist energy per unit length of the tube as U twist = k twist β , with β denoting the twist per unit length of the tube ( β = πατ ) and k twist denoting the torsionalstiffness. From the figure, we use the straight line fits near zero twist to estimate the valueof k twist as . eV nm and . eV nm, for the relaxed and unrelaxed cases respectively.It is also evident from the figure that non-linear effects start to play a noticeable role inthis nanotube at around degrees of twist per nanometer. These simulations serve as anexample of determining constitutive parameters directly from quantum mechanics usingHelical DFT. β ( degrees/nm ) T w i s t i n g e n e r g y p e r un i t l e n g t h ( e V / n m ) w/ atomic relaxtionStraight line fit near β = 0 w/o atomic relaxation Figure 7: Twist energy of a zigzag phosphorene nanotube as obtained via Helical DFT.
Since Helical DFT is an ab initio simulation tool, it can give us insights into the elec-tronic, optical and transport properties of materials under study. In particular, it may Due to the relatively large nanotube radius, the atoms in the untwisted nanotube are in an environ-ment similar to that in the phosphorene sheet. Since the phosphorene sheet itself had been relaxed well,the forces on the atoms in the nanotube were relatively small to begin with, and structural relaxationwas completed in just a few steps.
41e used to shed light into the coupling of mechanical strains with these properties. Toillustrate these points, we consider again the case of the zigzag phosphorene nanotubediscussed above (nanotube radius . nanometer, cyclic group order N = 64 .). Figure 8displays examples of helical band structure diagrams , which like their periodic counter-parts, can be used to illustrate the electronic levels in the system. A unique feature ofthese diagrams however, is that they can be used to display the variation of Kohn-Shameigenvalues with respect to η and ν . This makes them significantly easier to interpretthan traditional periodic band diagrams for quasi-one-dimesional (nanotube-like) struc-tures, which have only one index (i.e., the wave vector k z in the axial direction) labelingthe Kohn-Sham states. We anticipate that helical band diagrams like the ones shownare likely to emerge as powerful tools in understanding instabilities in optical and elec-tronic materials. The helical band diagrams described above can be used to compute the − . − . . . − . − . − . λ F η E n e r g y ( H a ) (a) Helical band diagram in η , along ν = 0 . − . − . − . λ F ν E n e r g y ( H a ) (b) Helical band diagram in ν , along η = 0 .Figure 8: Helical band structure diagrams for an untwisted zigzag phosphorene nanotube (radius = 3 . nanometers). λ F denotes the Fermi level (shifted to ).. size of the band gap of the system and infer whether it is conducting, semiconductingor insulating. We plot in Figure 9 the variation in the band gap of the above nan-otube, as it undergoes torsional deformations. The variation with and without atomicrelaxation effects are both displayed. From the figure, we observe that the nanotube hasa semiconducting behavior overall, and can be made to go from an insulating state atno twist (direct band gap of about . eV), to a practically conducting one, once thetwist reaches about degrees per nanometer. This is a rather significant change in theelectronic properties of the material, although in a mechanical sense, its deviation fromsimple linear elastic behavior (Figure 7) is still fairly modest at this level of twist. Similar band diagrams have been considered in [167, 168] in the context of phonon calculations ofcarbon nanotubes. For an extended system, such as a helical structure, the band gap is defined as the difference betweenthe conduction band minimum (CBM) and the valence band maximum (VBM). Within Helical DFT,this can be computed as the difference between the smallest eigenvalue above the Fermi level and thelargest eigenvalue below the Fermi level, as η is varied in I and ν in , , , . . . N − . Although it is well known that LDA is often unable to predict quantitatively accurate band gaps,the qualitative trends observed here are likely to be representative of actual physical behavior [169–173]in these nanotube systems. In any case, the formulation presented here does not have issues with regardto the use of more quantitatively accurate hybrid exchange correlation functionals [174, 175] whoseimplementation within our framework is a subject worthy of further investigation. . . . . β ( degrees/nm ) B a nd g a p ( e V ) w/ atomic relaxtionw/o atomic relaxation Figure 9: Effect of torsional deformation on the band gap of a zigzag phosphorene nanotube (radius = 3 . nanometers). To further illustrate the above electronic transition, we compute the (electronic) den-sity of states of the nanotube without and with twist ( β ≈ degrees/nanometer). Fol-lowing [43], we write this at a given energy level E and an electronic temperature of T e as: ℵ T e ( E ) = 2 (cid:90) I N N − (cid:88) ν =0 (cid:18) ∞ (cid:88) j =1 f (cid:48) T e (cid:0) E − λ j ( η, ν ) (cid:1)(cid:19) dη, (124)and evaluate it along a fine mesh of values of E in the range [ − , . The results areshown in Figure 10. It is apparent from the figure that the electronic states in the systemundergo significant change as the nanotube is subjected to twisting. In particular, thenumber of states at or near the Fermi level λ F , increases to values well above , indicatingonset of metallic behavior when the tube is twisted. − . − . − . λ F Energy (Ha) D e n s i t y o f s t a t e s ℵ T e ( · ) ( H a − ) (a) Density of states of untwisted nanotube. − . − . − . λ F ℵ T e ( λ F ) Energy (Ha) D e n s i t y o f s t a t e s ℵ T e ( · ) ( H a − ) w/ atomic relaxtionw/o atomic relaxation (b) Density of states of twisted nanotube.Figure 10: Electronic density of states plot for a zigzag phosphorene nanotube (radius = 3 . nanometers),without and with twist. λ F denotes the Fermi level (shifted to ). .5. Electronic properties: Behavior of chiral tubes under axial strains Finally, we use Helical DFT to study a chiral phosphorene nanotube. We choose atube with N = 64 and τ = 8 . Bohr as before, and also set α = 0 . . We relaxthe positions of the atoms within the simulation cell and use the resulting nanotube(observed to have an indirect band gap of about . eV) as the starting structure forsubsequent simulations. We subject this tube to (both tensile and compressive) axialstrains by varying τ , and relax the atomic positions in each case. Figure 11 shows thevariation in the band gap of this tube at different values of axial strain. The effect ofnot including atomic relaxation after the tube is subjected to the strains is also shown.From the figure, we see that in the range of strains considered, compressive strains tend − − . . . . . Percentage strain (%) B a nd g a p ( e V ) w/ atomic relaxtionw/o atomic relaxation Figure 11: Effect of axial strain on the band gap of a chiral phosphorene nanotube ( N = 64 , α = 0 . ). to reduce the band gap, while tensile strains appear to increase it. The tube appears totransition into a metallic state at about compressive strain. Based on the nature ofthe plot and motivated by earlier work on chiral carbon nanotubes [176, 177], we fitteda one-term Fourier series to this band gap data, and found that this produces a highquality fit for both the relaxed and unrelaxed cases. This suggests that like the case ofcarbon nanotubes, it might be possible to build (approximate) tight-binding type modelsof the band gap behavior [178] for phosphorene nanotubes. This is a topic that warrantsfurther investigation. Notably, the parameters in such models can be provided throughhigh-quality ab initio simulations based on Helical DFT.The above ab initio simulations of the phosphorene nanotubes (both zigzag and chiral)suggests that these materials have highly adjustable electronic states. Therefore, theymight find future applications as nanomaterials with tunable electronic/optical/transportproperties. The simulations also highlight the possibility of using inhomogeneous strainmodes for altering these properties, instead of homogeneous strain modes which are con-sidered in the strain-engineering literature [179–181] more commonly. Further investiga-tions of this material and others, along these lines, is the scope of future work.Finally, we find it worthwhile to point out that the phosphorene nanotube simulationsdescribed above (both electronic and mechanical) would be very challenging, or well-nighimpossible to carry out using conventional first principles techniques, even with the aidof massively parallel high-performance computing resources. In contrast, our MATLABimplementation of Helical DFT often allows such simulations to be carried out within44 few hours of simulation wall time on a desktop workstation or a single node of asupercomputing cluster, thus serving to reinforce the novelty and practical utility of theapproach.
5. Conclusions and Future Directions
In summary, we have presented Helical DFT — a novel, systematic first principles sim-ulation framework for systems with helical symmetries. We have presented a derivationof the equations of Kohn-Sham theory, as they apply to the case of helical structures. Ourderivation is systematic, self-contained and for the most part, mathematically rigorous.We have then solved these governing equations numerically by using a finite differencemethod based on helical coordinates. Using this working realization of the proposedapproach, we have carried out simulations involving phosphorene nanotubes, extractedtheir mechanical behavior ab initio , and identified changes in the electronic properties ofthis material as it undergoes twisting.Having laid out this foundational work on Helical DFT, we now discuss a number ofavenues for future investigation: • Development of an efficient spectral solution scheme: While the current finite dif-ference based implementation of Helical DFT allows us to investigate a number ofmaterials systems of interest, it also suffers a number of computational deficiencies(Section 3). Following our earlier work on the ab initio simulations of cluster sys-tems with arbitrary point group symmetries, we have already formulated, and arecurrently in the process of implementing a spectral scheme for solving the governingequations [8, 45]. This is expected to completely resolve the issues with the finitedifference formulation and provide a more efficient numerical implementation, thusopening the door to the simulation of more complex helical structures. • Mechanistic simulations of quasi-one-dimensional systems: The materials scienceliterature is rich with examples of quasi-one-dimensional structures that have beendiscovered and/or synthesized through experimental means. The simulation toolsdeveloped here can be used to characterize these materials computationally. Inparticular, the use of helical symmetry adapted ab initio molecular dynamics [182]can be used to study the mechanical behavior of these materials under axial andtorsional strains, as well as their instabilities and defects [48, 183]. • Helical Wannier states and the helical Berry phase: The helical Bloch Floquettransform U introduced in eqs. 37 and C.2 allows us to consider the notion ofhelical helical Wannier states. These, like their periodic counterparts, are definedthrough the action of U − on the helical Bloch states [82], and can be expectedto be exponentially localized for insulating systems [184]. Additionally, like thecase of periodic systems, a geometric phase (i.e., the Berry phase [185]) associatedwith the helical Bloch phase factor may be defined. These observations are likely tospur the development of novel computational analysis methods [186, 187] for helicalmaterials, as well as the discovery of novel topological materials [188, 189]. • Search for exotic materials properties and study of multi-physics coupling: SinceHelical DFT is a first principles simulation technique, it allows investigation of45he effect of torsional deformations on a material’s optical, electronic, magneticand transport properties. Like the case of the phospherene nanotube consideredin this work, it is possible that torsional or axial deformations in certain helicalstructures might induce a significant redistribution of electronic states in the ma-terial, leading to the appearance of interesting collective properties. In this regard,the investigation of the nanoscale flexoelectric effect [190–195] in helical systemswould be of particular interest, since owing to the quasi-one-dimensional nature ofthese materials, as well as the appearance of strain gradients in connection withtorsional deformations, a significant polarization may appear along the axis of ananostructure when twisted, thus leading to a strong flexoelectric effect. • Search for new phases of matter and coherent phase transformations: Finally, asremarked in [67, 68], the first principles techniques developed here might be in-strumental in the discovery of novel phases of matter. Additionally, the discoveryof transformations between such phases [67, 196] — particularly, coherent ones[197, 198] — are likely to lead to new classes of active materials. The methodsdeveloped here are likely to be very useful in the characterization of energy barriersof such transformations and help in their design.—
Appendix A. Verifying that the operator h aux η in Eq. 24 is symmetric The boundary conditions associated with the operator h aux η introduced in Eq. 24 arequite non-conventional. Here, we work through the steps of verifying that the operatoris symmetric when augmented with these boundary conditions. For this, we considersmooth functions f , f over D (the interior of the fundamental domain D G ), obeying thesame boundary conditions as φ in the proof of Theorem 2.3, i.e., f , ( x ) = f , (Υ h ◦ x ) , ∇ f , ( x ) = R − πα ∇ f , (Υ h ◦ x ) for x ∈ ∂ D z =0 and f , ( x ) = 0 for x ∈ ∂ D r = R . Then wehave: (cid:104) h aux η f , f (cid:105) L ( D ) = − (cid:90) D ∆ f f d x − i πητ (cid:90) D ∂f ∂x f d x − π η τ (cid:90) D f f d x + (cid:90) D V f f d x . (A.1)On using integration by parts [89], the first term on the right-hand side becomes: − (cid:90) D f ∆ f d x = 12 (cid:18) (cid:90) D ∇ f · ∇ f d x − (cid:90) ∂ D f ∇ f · d s (cid:19) , (A.2)where d s is the oriented surface measure. The boundary terms on the right-hand side ofeq. A.2 can be split as: (cid:90) ∂ D f ∇ f · d s = (cid:90) ∂ D r = R f ∇ f · d s + (cid:90) ∂ D z = τ f ∇ f · d s + (cid:90) ∂ D z =0 f ∇ f · d s . (A.3)The first term on the right-hand side of eq. A.3 goes to zero due to the Dirichlet bound-ary condition f ( x ) = 0 for x ∈ ∂ D r = R . To evaluate the second and the third terms,46ollectively denoted as I henceforth, we write the surface measures in terms of the localunit normals to get: I = (cid:90) ∂ D z = τ f ∇ f · ( e d s ) + (cid:90) ∂ D z =0 f ∇ f · ( − e d s )= (cid:90) ∂ D z = τ f ∇ f · ( e d s ) − (cid:90) ∂ D z =0 f ∇ f · ( e d s ) . (A.4)We now express y ∈ D z = τ as Υ h ◦ x with x ∈ D z =0 , so that the first integral on theright-hand side of eq. A.4 above can be rewritten as (with the dependence on x and y shown explicitly): (cid:90) y ∈ ∂ D z = τ f ( y ) ∇ f ( y ) · (cid:0) e d s y (cid:1) = (cid:90) x ∈ ∂ D z =0 f (Υ h ◦ x ) ∇ f (Υ h ◦ x ) · (cid:0) e d s Υ h ◦ x (cid:1) . (A.5)Here, the notation d s y is used to denote the surface measure centered at the point y , andsimilarly d s Υ h ◦ x denotes the surface measure centered at the point Υ h ◦ x . Since Υ h is anisometry, d s Υ h ◦ x has the same magnitude as d s x . Furthermore, the boundary conditionsimply that f (Υ h ◦ x ) = f ( x ) , and: ∇ f (Υ h ◦ x ) · e = (cid:0) R πα R − πα ∇ f (Υ h ◦ x ) (cid:1) · e = (cid:0) R − πα ∇ f (Υ h ◦ x ) (cid:1) · (cid:0) R T πα e (cid:1) = ∇ f ( x ) · e , (A.6)since R πα (and hence R T πα ) has axis e . Combining the above results, we arrive at: I = (cid:90) ∂ D z =0 f ( x ) ∇ f ( x ) · (cid:0) e d s x (cid:1) − (cid:90) ∂ D z =0 f ( x ) ∇ f ( x ) · (cid:0) e d s x (cid:1) = 0 , (A.7)and therefore, eq. A.2 reduces to: − (cid:90) D f ∆ f d x = 12 (cid:18) (cid:90) D ∇ f · ∇ f d x (cid:19) . (A.8)The second term on the right-hand side of eq. A.1 is: − i πητ (cid:20) (cid:90) D ∂f ∂x f d x (cid:21) = − i πητ (cid:20) (cid:90) D (cid:18) ∂f f ∂x − f ∂f ∂x (cid:19) d x (cid:21) . (A.9)We now consider the first integral on the right-hand side of eq. A.9, rewrite the integrandin divergence form, and use the Divergence Theorem [89] to get: (cid:90) D ∂f f ∂x d x = (cid:90) D ∇ · f f d x = (cid:90) ∂ D f f · d s , (A.10)with d s denoting the oriented surface measure (as earlier). The surface integral can besplit as: (cid:90) ∂ D f f · d s = (cid:90) ∂ D r = R f f · d s + (cid:90) ∂ D z = τ f f · d s + (cid:90) ∂ D z =0 f f · d s , (A.11)47rom which the first term on the right-hand side vanishes due to Dirichlet boundaryconditions on ∂ D r = R . For the second and the third terms, denoted collectively as I henceforth, we write the oriented surface measures in terms of the local unit normals toget: I = (cid:90) ∂ D z = τ f f · e d s − (cid:90) ∂ D z =0 f f · e d s = (cid:90) ∂ D z = τ f f d s − (cid:90) ∂ D z =0 f f d s (A.12)As earlier, we express y ∈ D z = τ as Υ h ◦ x with x ∈ D z =0 , so that the first integral on theright-hand side of eq. A.12 above can be rewritten by use of the boundary conditions as(the dependence on x and y has been shown explicitly): (cid:90) y ∈ ∂ D z = τ f ( y ) f ( y ) d s y = (cid:90) x ∈ ∂ D z =0 f (Υ h ◦ x ) f (Υ h ◦ y ) d s Υ h ◦ x = (cid:90) x ∈ ∂ D z =0 f ( x ) f ( x ) d s x . (A.13)It follows that, I = (cid:90) ∂ D z =0 f f d s − (cid:90) ∂ D z =0 f f d s = 0 , (A.14)and therefore, the second term on the right-hand side of eq. A.1 is: − i πητ (cid:20) (cid:90) D ∂f ∂x f d x (cid:21) = i πητ (cid:20) (cid:90) D f ∂f ∂x d x (cid:21) . (A.15)Combining eqs. A.1, A.8 and A.15, we get: (cid:104) h aux η f , f (cid:105) L ( D ) =12 (cid:90) D ∇ f · ∇ f d x + i πητ (cid:90) D f ∂f ∂x d x − π η τ (cid:90) D f f d x + (cid:90) D V f f d x . (A.16)On the other hand, we can express (cid:104) f , h aux η f (cid:105) L ( D ) as: (cid:90) D f ( −
12 ∆ f ) d x + (cid:90) D f (cid:18) − i πητ ∂f ∂x (cid:19) d x + (cid:90) D f (cid:18) − π η τ f (cid:19) d x + (cid:90) D f (cid:0) V f (cid:1) d x = − (cid:90) D f ∆ f d x + i πητ (cid:90) D f ∂f ∂x d x − π η τ (cid:90) D f f d x + (cid:90) D V f f d x . (A.17)Integrating by parts the first term on the right-hand side of eq. A.17, we have: − (cid:90) D f ∆ f d x = 12 (cid:18) (cid:90) D ∇ f · ∇ f d x − (cid:90) ∂ D f ∇ f · d s (cid:19) . (A.18)48he second term on the right-hand side of eq. A.18 can be sent to zero by applicationof the boundary conditions (using a procedure similar to the one outlined in eqs. A.3 -A.7). This leaves us with: (cid:104) f , h aux η f (cid:105) L ( D ) = 12 (cid:90) D ∇ f · ∇ f + i πητ (cid:90) D f ∂f ∂x d x − π η τ (cid:90) D f f d x + (cid:90) D V f f d x = (cid:104) h aux η f , f (cid:105) L ( D ) , (A.19)which implies that the operator h aux η is symmetric on smooth functions obeying the bound-ary conditions outlined above. Since such functions are dense in the domain of h aux η (theboundary conditions being interpreted in the trace sense in that case), the symmetry ofthe operator follows. Appendix B. Direct integral decomposition of the single electron Hamilto-nian
The formalism of direct integrals [85, 104] generalizes the idea of direct sums in Hilbertspaces and it allows us to make the idea of diagonalizing or block-diagonalizing an un-bounded self-adjoint operator mathematically precise. Here we provide a brief summaryof some of the key ideas associated with direct integrals. We then demonstrate how thehelical Bloch-Floquet transform can be employed to effectively “block-diagonalize” thesingle electron Hamiltonian in the sense of direct integrals. This can be viewed as anatural extension of the block-diagonal decomposition of the single-electron Hamiltonianin the sense of direct sums, that applies when structures associated with finite symmetrygroups are considered [8, 39].If H (cid:48) is a (separable) Hilbert space and ( M, µ ) is a ( σ -finite) measure space, thenthe Hilbert space H = L ( M, dµ ; H (cid:48) ) of H (cid:48) valued functions which are square integrable(against the measure µ ), is defined to be a (constant fiber) direct integral , and wedenote this relationship as: H = (cid:90) ⊕ M H (cid:48) dµ . (B.1)Vector addition and scalar multiplication are defined pointwise in this space, i.e., for s ∈ M , f , f ∈ H and z ∈ C , we have (cid:0) f + f (cid:1) ( s ) = f ( s ) + f ( s ) and (cid:0) z f (cid:1) ( s ) = z f ( s ) ,while the inner product is defined as: (cid:104) f , f (cid:105) H = (cid:90) M (cid:104) f ( s ) , f ( s ) (cid:105) H (cid:48) dµ . (B.2)The above definitions allows us to decompose operators on H in terms of operators on H (cid:48) in a particular sense. Let L ( H (cid:48) ) denote the set of the bounded linear operators on This definition subsumes the notion of direct sums: if µ is a sum of point measures at a fi-nite set of points s , s , . . . , s K , then any f ∈ L ( M, dµ ; H (cid:48) ) is determined by the set of K values { f ( s ) , f ( s ) , . . . , f ( s K ) } . Thus, H = L ( M, dµ ; H (cid:48) ) is isomorphic to the direct sum K (cid:77) i =1 H (cid:48) . H (cid:48) . A bounded linear operator A on the space H (cid:18) = (cid:90) ⊕ M H (cid:48) dµ (cid:19) is said tobe decomposed in the sense of a direct integral decomposition if there exists a function A ( · ) in L ∞ (cid:0) M, dµ ; L ( H (cid:48) ) (cid:1) such that for all f ∈ H , the relationship ( A f )( s ) = A ( s ) f ( s ) holds. We then call A decomposable , we refer to the operators A ( s ) as the fibers of A ,and we denote this relationship as: A = (cid:90) ⊕ M A ( s ) dµ . (B.3)Conversely, it is also possible to “build” operators on the space H , starting fromoperators on the space H (cid:48) . Specifically, given a measurable function A ( · ) from M tothe set of (bounded or unbounded) self-adjoint operators the on H (cid:48) , we may define anoperator A on H with domain:Dom. ( A ) = (cid:26) f ∈ H : f ( s ) ∈ D ( A ( s )) a.e. ; (cid:90) M (cid:107) A ( s ) f ( s ) (cid:107) H (cid:48) dµ < ∞ (cid:27) , (B.4)as: ( A f )( s ) = A ( s ) f ( s ) . (B.5)As before, we will use the notation: A = (cid:90) ⊕ M A ( s ) dµ , (B.6)to denote the above relationship between the operator A and its fibers A ( s ) .With these definitions in place, we now apply the above apparatus to carry out asuitable decomposition of the single-electron Hamiltonian associated with a helical struc-ture. We will consider the case of a structure associated with a helical group G that isgenerated by a single element Υ h . As discussed in 2.2, the single-electron Hamiltonianin this case is the operator H = − ∆ + V ( x ) over the space (cid:101) H = L ( C ) . The relevantboundary condition is ψ ( x ) = 0 for x ∈ ∂ C and the potential V ( x ) is group invariant,i.e., V ( x ) = V ( (cid:101) Υ ◦ x ) for every (cid:101) Υ ∈ G .We may choose M = [ − , ) = I , H (cid:48) = L ( D ) , µ as the Lebesgue measure on I (denoted as dη henceforth). Then, H = (cid:90) ⊕ I H (cid:48) dη . (B.7)is the space L ( D × I ) . The helical Bloch-Floquet transform U introduced in eq. 37 allowsus to map functions and operators in between the spaces (cid:101) H = L ( C ) and H = L ( D × I ) .Specifically, for each η ∈ I , let H η be the restriction of H to D along with the boundaryconditions ψ (Υ h ◦ x ) = e − i πη ψ ( x ) , R − πα ∇ ψ (Υ h ◦ x ) = e − i πη ∇ ψ ( x ) for x ∈ ∂ D z =0 A function A ( · ) from M to L ( H (cid:48) ) is called measurable if for every f , f ∈ H (cid:48) , the map s (cid:55)→(cid:104) f , A ( s ) f (cid:105) H (cid:48) is measurable [85]. A function A ( · ) from M to the set of (bounded or unbounded) self-adjoint operators on H (cid:48) is calledmeasurable if the function ( A ( · ) + i ) − is measurable [85]. ψ ( x ) = 0 for x ∈ ∂ D r = R . Then, the operator H can be decomposed into the set ofoperators (cid:8) H η (cid:9) η ∈ I (referred to as the fibers of H ) using the direct integral decomposition,in the sense: U H U − = (cid:90) ⊕ I H η dη , (B.8)To demonstrate this result, we will establish direct integral decompositions of the Lapla-cian and the potential operator individually as operators on L ( C ) , i.e., we will showsuccessively that: U (cid:0) − ∆ (cid:1) U − = (cid:90) ⊕ I ( − ∆) η dη , (B.9)and, U V U − = (cid:90) ⊕ I V η dη . (B.10)The sought result will then follow from Theorem XIII.85g of [85].First, to establish eq. B.9, i.e., to perform the direct integral decomposition of theLaplacian on L ( C ) , subject to zero Dirichlet boundary condition for x ∈ C , we definethe fibers ( − ∆) η as the restriction of − ∆ to D along with the boundary conditions f (Υ h ◦ x ) = e − i πη f ( x ) , R − πα ∇ f (Υ h ◦ x ) = e − i πη ∇ f ( x ) for x ∈ ∂ D z =0 and, f ( x ) = 0 for x ∈ ∂ D r = R . Let B be the operator on the right-hand side of eq. B.9 acting on thespace H = L ( D × I ) , and let f be a Schwartz class function that obeys f ( x ) = 0 for x ∈ C . Then U f is given as: ( U f )( x , η ) = (cid:88) m ∈ Z f (Υ m h ◦ x ) e i πmη , (B.11)and it is a smooth function in x over D that obeys ( U f )( x , η ) = 0 for x ∈ ∂ D r = R forany η ∈ I . Furthermore, since f is a Schwartz class function, so are its derivatives, andtherefore we may take derivatives of the above expression to get that: ∂ n ( U f )( x , η ) ∂x ni = U (cid:18) ∂ n f∂x ni (cid:19) ( x , η ) , (B.12)for every η ∈ I , any n ∈ N and i = 1 , , . From this, and using the series expressionfor U f , it easily follows that the required boundary conditions f (Υ h ◦ x ) = e − i πη f ( x ) , R − πα ∇ f (Υ h ◦ x ) = e − i πη ∇ f ( x ) for x ∈ ∂ D z =0 are satisfied. Thus for every η ∈ I , it holdsthat U f ( · ; η ) ∈ Dom. (cid:0) ( − ∆) η (cid:1) and that − ∆ η ( U f ( · ; η )) = U ( − ∆)( · ; η ) . Thus, it followsthat U f ∈ Dom. ( B ) and the relationship B ( U f ) = U ( − ∆) holds, which establishes B.9.Next, to establish B.10, we set the fibered operator V η on L ( D ) as: ( V η f )( x ) = V ( x ) f ( x ) . (B.13)Then for any Schwartz class function f over C , we have that: U ( V f )( x ; η ) = (cid:88) m ∈ Z V (Υ m h ◦ x ) f (Υ m h ◦ x ) e i πmη , (B.14)51hich, using the invariance of the potential under G , gives: U ( V f )( x ; η ) = V ( x ) (cid:88) m ∈ Z f (Υ m h ◦ x ) e i πmη = V ( x ) U ( x ; η ) = V η ( U f )( x , η ) . (B.15)Thus B.10 is established.With B.8 established, some conclusions can be immediately drawn regarding the struc-ture of the spectrum of H . Using Theorem XIII.85 in [85] for example, it follows thatthe collection of eigenvalues of the operators H η together form the spectrum of H , i.e., Λ = spec . ( H ) . Furthermore, Theorem XIII.86 of [85] implies that H has a purely contin-uous spectrum.Finally, a result made use of in Section 2.2.2 is that if A is an operator that is invariantunder G (i.e., it commutes with the unitary operators in the set T = (cid:8) T (cid:101) Υ : T (cid:101) Υ f ( x ) = f ( (cid:101) Υ − ◦ x ) (cid:9) (cid:101) Υ ∈G ) and it is is locally trace-class, we may assign meaning to the trace perunit fundamental domain as: Tr. [ A ] = (cid:90) I Tr. [ A η ] dη . (B.16)Discussion and rigorous proofs of this result appear in [199, 200] for the case of periodicsymmetries. To see why this result also holds true for the case of a helical symmetry group,we denote I D as the indicator function of the fundamental domain. Let { f j : j ∈ N } bean orthonormal basis of L ( C ) . Due to the fact that A is invariant under the group G ,we may write: U A U − = (cid:90) ⊕ I A η dη , (B.17)Now, the trace per unit fundamental domain is:Tr. [ A ] = ∞ (cid:88) j =1 (cid:104) I D A I D f j , f j (cid:105) L ( C ) = ∞ (cid:88) j =1 (cid:104)A I D f j , I D f j (cid:105) L ( C ) (B.18)Using eq. B.17, this becomes:Tr. [ A ] = ∞ (cid:88) j =1 (cid:28) U − (cid:18) (cid:90) ⊕ I A η dη (cid:19) U I D f j , I D f j (cid:29) L ( C ) = ∞ (cid:88) j =1 (cid:28)(cid:18) (cid:90) ⊕ I A η dη (cid:19) U I D f j , U I D f j (cid:29) L ( C ) , (B.19)which, upon using B.5 and B.2, can be written as:Tr. [ A ] = ∞ (cid:88) j =1 (cid:90) I (cid:104)A η U I D f j , U I D f j (cid:105) L ( D ) dη (B.20)Now, recognizing that the indicator function I D acts as a projection operator on L ( C ) , wesee that U I D f j are simply basis functions of L ( D ) for every η ∈ I . Denoting ˜ f j ( x ; η ) =( U I D f j )( x ; η ) and exchanging the summation and the integral, the above equation canbe re-written as:Tr. [ A ] = (cid:90) I ∞ (cid:88) j =1 (cid:104)A η ˜ f j ( · ; η ) , ˜ f j ( · ; η ) (cid:105) L ( D ) dη = (cid:90) I Tr. [ A η ] dη , (B.21)as claimed. 52 ppendix C. Helical structure associated with a group generated by twogroup elements: Expressions for important physical quantitiesand governing equations For a helical structure associated with a group G that is generated by two groupelements (i.e., a screw transformation of the form ( R πmα | mτ e ) and a pure rotationof the form ( R n Θ | ) with Θ = π N , and both rotations R πmα , R n Θ having e as thecommon axis) the form of the governing equations as well as the expressions for thevarious quantities of interest can be easily deduced by using the discussion in Section2.2 as a starting point, and using the following rules of substitution. The helical Blochphase factors are of the form e − i π ( mη + nν N ) with n, ν ∈ { , , . . . , N − } and η ∈ I ,instead of being e − i πmη only. Consequently, helical Bloch states and the helical bandsare labeled as ψ j ( x ; η, ν ) and λ j ( η, ν ) , respectively. Integrals over η ∈ I = [ − , ) need tobe replaced with integrals over η ∈ I and summations over ν of the form N N − (cid:88) ν =0 . Takinginto account the effect of all the symmetry operations of the group (i.e., computing thegroup orbit of a point for example) amounts to summing over all group elements of theform Υ m h • Υ n c = ( R πmα + n Θ | mτ e ) with m ∈ Z and n ∈ { , , . . . , N − } . Also, in theexpression for the forces (eq. 107), the rotation matrix ( R πmα ) − needs to be replacedwith ( R πmα + n Θ ) − . Finally, spatial integrals over the fundamental domain D G (or itsinterior D ) have to be replaced by integrals over the fundamental domain D G (or itsinterior (cid:101) D ). The reciprocal space (or more specifically, the Brillouin zone of the reciprocalspace) for the symmetry group G will be denoted by the set B = I × { , , , . . . , N − } .The mathematical reasons behind the above substitution rules are as follows: Thestructure of the helical group generated by two elements is such that it can be expressedas the direct product of the helical group generated by a single element and a cyclic group.Therefore, the Bloch theorem for a structure associated with such a symmetry group canbe established by first using Theorem 2.3, and then using an appropriate version ofTheorem 2.6 from [39]. The characters (i.e., one-dimensional complex irreducible repre-sentations) of the helical group generated by two elements can be obtained by multiplyingout the characters of the helical group generated by a single element and the cyclic group.This leads to the helical Bloch states in this case to obey the condition: ψ j (cid:0) (Υ m h • Υ n c ) ◦ x ; η, ν (cid:1) = e − i π ( mη + nν N ) ψ j ( x ; η, ν ) . (C.1)The completeness of these states follows by a combination of Theorem 2.6 and the prop-erties of the Peter-Weyl projectors [39–41] for the cyclic group. This completeness resultpaves the way for a suitable Helical Bloch-Floquet transform U : L ( C ) → L ( (cid:101) D × B ) : ( U f )( x , η, ν ) = (cid:88) m ∈ Z (cid:18) N N − (cid:88) n =0 f (cid:0) (Υ m h • Υ n c ) ◦ x (cid:1) e i π ( mη + nν N ) (cid:19) , (C.2)as well as a direct integral representation of the single electron Hamiltonian. From these,the key quantities of interest and the governing equations may be systematically deduced Theorem 2.6 from [39] was established for the single electron Hamiltonian operator. However, it canbe easily extended to any other linear elliptic self-adjoint operator that commutes the cyclic symmetriesof the system, e.g. h aux η defined in eq. 24, with V ( x ) for that operator invariant under G . G . For the sakeof definiteness, we assume that the helical structure, is embedded in the infinite cylinder C (eq. 10), and as in eq. 11, we let D G = (cid:8) ( r, ϑ, z ) : 0 ≤ r < R, ≤ ϑ < π N , ≤ z < τ (cid:9) denote a cylindrical sector that acts as the fundamental domain of the helical structure.Let (cid:101) D denote the interior of this fundamental domain. Furthermore, let the positions ofthe atoms within the fundamental domain be P G = (cid:8) x k (cid:9) M G k =1 and let (cid:8) b k ( x ; x k ) (cid:9) M G k =1 de-note the corresponding nuclear pseudocharges. The helical Bloch states presented abovein eq. C.1 allow the single electron problem for the entire helical structure to be reducedto the fundamental domain D G (or its interior (cid:101) D ).Let g j ( η, ν ) = f T e (cid:0) λ j ( η, ν ) (cid:1) denote the electronic occupation numbers at electronictemperature T e . The electron density for x ∈ (cid:101) D can be expressed as: ρ ( x ) = 2 (cid:90) I N N − (cid:88) ν =0 ∞ (cid:88) j =1 g j ( η, ν ) | ψ j ( x ; η, ν ) | dη . (C.3)The electron density needs to obey the constraint of integrating to a fixed number ofelectrons in the fundamental domain, i.e., (cid:90) (cid:101) D ρ ( x ) d x = N e = ⇒ (cid:90) I N N − (cid:88) ν =0 ∞ (cid:88) j =1 g j ( η, ν ) = N e . (C.4)The electronic free energy per unit fundamental domain can be expressed as: F (Λ , Ψ , P G , (cid:101) D, G ) = T s (Λ , Ψ , (cid:101) D, G ) + E xc ( ρ, (cid:101) D ) + K (Λ , Ψ , P G , (cid:101) D, G ) (C.5) + E el ( ρ, P G , (cid:101) D, G ) − T e S (Λ) . (C.6)The terms on the right-hand side of the above equation are as follows. The first term isthe kinetic energy of the electrons per unit fundamental domain and is expressible as: T s (Λ , Ψ , (cid:101) D , G ) = (cid:90) I N N − (cid:88) ν =0 (cid:18) ∞ (cid:88) j =1 g j ( η, ν ) (cid:10) ∆ ψ j ( · ; η, ν ) , ψ j ( · ; η, ν ) (cid:11) L ( (cid:101) D ) dη (cid:19) . (C.7)The second term is the exchange correlation energy per unit fundamental domain and isexpressible as: E xc ( ρ, (cid:101) D ) = (cid:90) (cid:101) D ε xc [ ρ ( x )] ρ ( x ) d x . (C.8)The third term is the nonlocal pseudopotential energy per unit fundamental domain andis expressible as: K (Λ , Ψ , P G , (cid:101) D , G )= 2 M G (cid:88) k =1 (cid:88) p ∈N k ∞ (cid:88) j =1 γ k,p (cid:90) I N N − (cid:88) ν =0 (cid:18) g j ( η, ν ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:10) ˆ χ k,p ( · ; η, ν ; x k ) , ψ j ( · ; η, ν ) (cid:11) L ( (cid:101) D ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:19) dη , (C.9)54ith the functions: ˆ χ k,p ( x ; η, ν ; x k ) = (cid:88) m ∈ Z N − (cid:88) n =0 χ k ; p (cid:0) x ; (Υ m h • Υ n c ) ◦ x k (cid:1) e − i π ( mη + nν N ) , (C.10)Here, as in eq. 76, χ k ; p ( · ; p k ) denote the atom centered nonlocal pseudopotential oper-ator projection functions. The fourth term on the right-hand side of eq. C.6 is the netelectrostatic interaction energy per unit fundamental domain, and can be written as: E el ( ρ, P G , (cid:101) D , G ) = max Φ (cid:26) − π (cid:90) (cid:101) D |∇ Φ | d x + (cid:90) (cid:101) D (cid:18) ρ ( x ) + b ( x , , P G , G ) (cid:19) Φ( x ) d x (cid:27) + E sc ( P G , G , (cid:101) D ) . (C.11)Here, b ( x , , P G , G ) is the net nuclear pseudocharge: b ( x , P G , G ) = (cid:88) m ∈ Z N − (cid:88) n =0 M G (cid:88) k =1 b k (cid:0) x ; (Υ m h • Υ n c ) ◦ x k (cid:1) , (C.12)and E sc ( P G , G , (cid:101) D ) represents self interaction and correction terms. The net electrostaticpotential Φ can be expressed as the Newtonian potential: Φ( x ) = (cid:90) C ρ ( y ) + b ( y , P G , G ) (cid:107) x − y (cid:107) R d y , (C.13)The last term on the right-hand side of eq. C.6 is the electronic entropy contribution: S (Λ)= − k B (cid:90) I N N − (cid:88) ν =0 (cid:20) ∞ (cid:88) j =1 g j ( η, ν ) log (cid:0) g j ( η, ν ) (cid:1) + (cid:0) − g j ( η, ν ) (cid:1) log (cid:0) − g j ( η, ν ) (cid:1)(cid:21) dη (C.14)The symmetry adapted Kohn-Sham equations for the system, as posed over the fun-damental domain are: H KS [Λ , Ψ , P G , (cid:101) D, G ] ψ j ( · ; η, ν ) = λ j ( η, ν ) ψ j ( · ; η, ν ) , (C.15)for j ∈ N , η ∈ I , ν ∈ { , , . . . , N − } . Here, H KS [Λ , Ψ , P G , (cid:101) D, G ] ≡ −
12 ∆ + V xc + Φ + (cid:98) V nl (cid:101) D , (C.16)is the Kohn-Sham operator, with its dependence on the helical Bloch sates, the helicalbands, etc., made explicit . As before, V xc represents the exchange correlation potential, As in the discussion in Section 2.2.2, we may change notation and express the free energy per unit cellin terms of the occupations as (cid:101) F ( G , Ψ , P G , (cid:101) D , G ) = F (Λ , Ψ , P G , (cid:101) D , G ) , instead of using the eigenvalues.With this change in notation, the problem of determination of the Kohn-Sham ground state for a givenhelical structure, can be expressed as: (cid:101) F ( P G , (cid:101) D , G ) = inf. Ψ , G (cid:101) F ( G , Ψ , P G , (cid:101) D , G ) , (C.17)subject to the condition that Ψ consists of helical Bloch states and eq. C.4. Similarly, the Kohn Shamoperator in eqs. C.15 and C.16 may be denoted as as H KS [ G , Ψ , P G , (cid:101) D, G ] . (cid:98) V nl (cid:101) D = M G (cid:88) k =1 (cid:88) p ∈N k γ k,p ˆ χ k,p ( · ; η, ν ; x k ) ˆ χ k,p ( · ; η, ν ; x k ) (C.18)The helical Bloch states obey (for i, j ∈ N , ): (cid:104) ψ i ( · ; η, ν ) , ψ j ( · ; η, ν ) (cid:105) L ( (cid:101) D ) = δ i,j . (C.19)The net electrostatic potential obeys: − ∆Φ = 4 π (cid:0) ρ + b ( · , P G , G ) (cid:1) , (C.20)and is invariant under G . Finally, the Harris-Foulkes energy functional is: F HF (Λ , Ψ , P G , (cid:101) D , G ) = 2 (cid:90) I N N − (cid:88) ν =0 ∞ (cid:88) j =1 λ j ( η ) g j ( η ) dη + E xc ( ρ, (cid:101) D ) − (cid:90) (cid:101) D V xc ( ρ ( x )) ρ ( x ) d x + 12 (cid:90) (cid:101) D (cid:18) b ( x , , P G , G ) − ρ ( x ) (cid:19) Φ( x ) d x + E sc ( P G , G , (cid:101) D ) − T e S (Λ) . (C.21)and the Hellman-Feynman forces on the atoms in the fundamental domain are: f k = (cid:88) m ∈ Z N − (cid:88) n =0 ( R πmα + n Θ ) − (cid:90) (cid:101) D ∇ b k (cid:0) x ; (Υ m h • Υ n c ) ◦ x k (cid:1) Φ( x ) d x − ∂E sc ( P G , G , (cid:101) D ) ∂ x k − ∞ (cid:88) j =1 (cid:32) (cid:90) I N N − (cid:88) ν =0 g j ( η, ν ) (cid:88) p ∈N k γ k ; p Re. (cid:40)(cid:20) (cid:90) (cid:101) D ˆ χ k,p ( x ; η, ν ; x k ) ψ j ( x ; η, ν ) d x (cid:21) × (cid:20) (cid:90) (cid:101) D ψ j ( x ; η, ν ) ∂ ˆ χ k,p ( x ; η, ν ; x k ) ∂ x k d x (cid:21)(cid:41)(cid:33) dη (C.22)A more comprehensive account of the above equations as well as extensive numericalsimulations involving them appears in a forthcoming contribution [70]. Appendix D. Cartesian gradient operator, Laplacian operator and integralsin helical coordinates
We derive expressions in helical coordinates for the Cartesian gradient operator (usefulin expressing forces on the atoms), the Laplacian operator (useful for expressing theSchrodinger operator), and integrals (useful for computing energies) in this Appendix.The calculations involved are straight forward applications of the chain rule and theyoriginally appear in [8]. We include these here for the sake of completeness.Let ξ ( x , x , x ) be a generic scalar quantity expressed in Cartesian coordinates. Inhelical coordinates, this can be expressed as: ξ ( x , x , x ) = ˆ ξ (cid:0) r ( x , x , x ) , θ ( x , x , x ) , θ ( x , x , x ) (cid:1) , (D.1)56e wish to evaluate the Cartesian gradient: ∇ ξ = ∂ξ∂x e + ∂ξ∂x e + ∂ξ∂x e , (D.2)and the Laplacian: ∆ ξ = ∂ ξ∂x + ∂ ξ∂x + ∂ ξ∂x , (D.3)in terms of the function ˆ ξ and the coordinates r, θ , θ . Let i, j = 1 , . . . , . We then haveby the chain rule: ∂ξ∂x i = ∂ ˆ ξ∂r ∂r∂x i + ∂ ˆ ξ∂θ ∂θ ∂x i + ∂ ˆ ξ∂θ ∂θ ∂x i , (D.4)and further, ∂ ξ∂x i ∂x j = ∂ ˆ ξ∂r ∂ r∂x i ∂x j + ∂ ˆ ξ∂θ ∂ θ ∂x i ∂x j + ∂ ˆ ξ∂θ ∂ θ ∂x i ∂x j + ∂ (cid:0) ∂ ˆ ξ∂r (cid:1) ∂x j ∂r∂x i + ∂ (cid:0) ∂ ˆ ξ∂θ (cid:1) ∂x j ∂θ ∂x i + ∂ (cid:0) ∂ ˆ ξ∂θ (cid:1) ∂x j ∂θ ∂x i . (D.5)The chain rule applied again to the last 3 terms gives us: ∂ (cid:0) ∂ ˆ ξ∂r (cid:1) ∂x j ∂r∂x i = ∂r∂x i (cid:18) ∂ ˆ ξ∂r ∂r∂x j + ∂ ˆ ξ∂r∂θ ∂θ ∂x j + ∂ ˆ ξ∂r∂θ ∂θ ∂x j (cid:19) .∂ (cid:0) ∂ ˆ ξ∂θ (cid:1) ∂x j ∂θ ∂x i = ∂θ ∂x i (cid:18) ∂ ˆ ξ∂θ ∂r ∂r∂x j + ∂ ˆ ξ∂θ ∂θ ∂x j + ∂ ˆ ξ∂θ ∂θ ∂θ ∂x j (cid:19) .∂ (cid:0) ∂ ˆ ξ∂θ (cid:1) ∂x j ∂θ ∂x i = ∂θ ∂x i (cid:18) ∂ ˆ ξ∂θ ∂r ∂r∂x j + ∂ ˆ ξ∂θ ∂θ ∂θ ∂x j + ∂ ˆ ξ∂θ ∂θ ∂x j (cid:19) . (D.6)For i = j , the above expressions can be combined to yield: ∂ ˆ ξ∂x i = ∂ ˆ ξ∂r (cid:18) ∂r∂x i (cid:19) + ∂ ˆ ξ∂θ (cid:18) ∂θ ∂x i (cid:19) + ∂ ˆ ξ∂θ (cid:18) ∂θ ∂x i (cid:19) + 2 (cid:18) ∂ ˆ ξ∂r∂θ ∂r∂x i ∂θ ∂x i + ∂ ˆ ξ∂θ ∂θ ∂θ ∂x i ∂θ ∂x i + ∂ ˆ ξ∂θ ∂r ∂θ ∂x i ∂r∂x i (cid:19) . (D.7)The helical coordinates and their first derivatives with respect to the Cartesian coordi-nates x , x , x are as follows: r ( x , x , x ) = (cid:113) x + x , θ ( x , x , x ) = x τ , θ ( x , x , x ) = 12 π arctan ( x x ) − α x τ ,∂r∂x = x (cid:112) x + x = x r , ∂r∂x = x (cid:112) x + x = x r , ∂r∂x = 0 ,∂θ ∂x = 0 , ∂θ ∂x = 0 , ∂θ ∂x = 1 τ ,∂θ ∂x = − π x x + x = − π x r , ∂θ ∂x = 12 π x x + x = 12 π x r , ∂θ ∂x = − ατ . (D.8)57sing these, we get: ∂ξ∂x = ∂ ˆ ξ∂r ∂r∂x + ∂ ˆ ξ∂θ ∂θ ∂x + ∂ ˆ ξ∂θ ∂θ ∂x = ˆ ξ r cos (cid:0) π ( αθ + θ ) (cid:1) − ˆ ξ θ sin (cid:0) π ( αθ + θ ) (cid:1) πr , (D.9) ∂ξ∂x = ∂ ˆ ξ∂r ∂r∂x + ∂ ˆ ξ∂θ ∂θ ∂x + ∂ ˆ ξ∂θ ∂θ ∂x = ˆ ξ r sin (cid:0) π ( αθ + θ ) (cid:1) + ˆ ξ θ cos (cid:0) π ( αθ + θ ) (cid:1) πr , (D.10) ∂ξ∂x = ∂ ˆ ξ∂r ∂r∂x + ∂ ˆ ξ∂θ ∂θ ∂x + ∂ ˆ ξ∂θ ∂θ ∂x = 1 τ (cid:0) ˆ ξ θ − α ˆ ξ θ (cid:1) , (D.11)which completes the calculation of the Cartesian gradient.We now compute the second derivatives of the helical coordinates with respect to theCartesian coordinates, but using eq. D.7, we restrict ourselves only to the quantities thatwould appear in the Laplacian: ∂ r∂x = 1 r − x r , ∂ θ ∂x = 0 , ∂ θ ∂x = 12 π x x r ,∂ r∂x = 1 r − x r , ∂ θ ∂x = 0 , ∂ θ ∂x = − π x x r ,∂ r∂x = 0 , ∂ θ ∂x = 0 , ∂ θ ∂x = 0 . (D.12)We are now ready to evaluate (D.5) through (D.7), (D.8) and (D.12): ∂ ξ∂x = ∂ ˆ ξ∂r ∂ r∂x + ∂ ˆ ξ∂θ ∂ θ ∂x + ∂ ˆ ξ∂θ ∂ θ ∂x + ∂ ˆ ξ∂r (cid:18) ∂r∂x (cid:19) + ∂ ˆ ξ∂θ (cid:18) ∂θ ∂x (cid:19) + ∂ ˆ ξ∂θ (cid:18) ∂θ ∂x (cid:19) + 2 (cid:18) ∂ ˆ ξ∂r∂θ ∂r∂x ∂θ ∂x + ∂ ˆ ξ∂θ ∂θ ∂θ ∂x ∂θ ∂x + ∂ ˆ ξ∂θ ∂r ∂θ ∂x ∂r∂x (cid:19) = ˆ ξ r (cid:18) r − x r (cid:19) + ˆ ξ θ (cid:18) x x πr (cid:19) + ˆ ξ rr x r + ˆ ξ θ θ (cid:18) x π r (cid:19) − ˆ ξ θ r (cid:18) x x πr (cid:19) . (D.13)Similarly, ∂ ξ∂x = ∂ ˆ ξ∂r ∂ r∂x + ∂ ˆ ξ∂θ ∂ θ ∂x + ∂ ˆ ξ∂θ ∂ θ ∂x + ∂ ˆ ξ∂r (cid:18) ∂r∂x (cid:19) + ∂ ˆ ξ∂θ (cid:18) ∂θ ∂x (cid:19) + ∂ ˆ ξ∂θ (cid:18) ∂θ ∂x (cid:19) + 2 (cid:18) ∂ ˆ ξ∂r∂θ ∂r∂x ∂θ ∂x + ∂ ˆ ξ∂θ ∂θ ∂θ ∂x ∂θ ∂x + ∂ ˆ ξ∂θ ∂r ∂θ ∂x ∂r∂x (cid:19) = ˆ ξ r (cid:18) r − x r (cid:19) − ˆ ξ θ (cid:18) x x πr (cid:19) + ˆ ξ rr x r + ˆ ξ θ θ x π r + ˆ ξ θ r x x πr . (D.14)58nd, ∂ ξ∂x = ∂ ˆ ξ∂r ∂ r∂x + ∂ ˆ ξ∂θ ∂ θ ∂x + ∂ ˆ ξ∂θ ∂ θ ∂x + ∂ ˆ ξ∂r ( ∂r∂x ) + ∂ ˆ ξ∂θ (cid:18) ∂θ ∂x (cid:19) + ∂ ˆ ξ∂θ (cid:18) ∂θ ∂x (cid:19) + 2 (cid:18) ∂ ˆ ξ∂r∂θ ∂r∂x ∂θ ∂x + ∂ ˆ ξ∂θ ∂θ ∂θ ∂x ∂θ ∂x + ∂ ˆ ξ∂θ ∂r ∂θ ∂x ∂r∂x (cid:19) = ˆ ξ θ θ τ + ˆ ξ θ θ α τ − ξ θ θ ατ . (D.15)So, we have: ∆ ξ = ∂ ξ∂x + ∂ ξ∂x + ∂ ξ∂x = ˆ ξ rr + 1 r ˆ ξ r + 1 τ ˆ ξ θ θ − ατ ˆ ξ θ θ + 14 π (cid:18) r + 4 π α τ (cid:19) ˆ ξ θ θ . (D.16)Finally, we compute the Jacobian determinant of the transformation to helical coor-dinates as: = Det. ∂x ∂r ∂x ∂θ ∂x ∂θ ∂x ∂r ∂x ∂θ ∂x ∂θ ∂x ∂r ∂x ∂θ ∂x ∂θ = Det. cos(2 π ( αθ + θ )) − παr sin(2 π ( αθ + θ )) − πr sin(2 π ( αθ + θ ))sin(2 π ( αθ + θ )) 2 παr cos(2 π ( αθ + θ )) 2 πr cos(2 π ( αθ + θ ))0 τ = 2 πτ r (D.17)Thus, the integral of a scalar function ξ ( x ) over the simulation cell Ω (expressed incylindrical coordinates as Ω = (cid:8) ( r, ϑ, z ) ∈ R : R in ≤ r ≤ R out , ≤ ϑ ≤ Θ , ≤ z ≤ τ (cid:9) )can be written as: (cid:90) Ω ξ ( x ) d x = (cid:90) (cid:90) (cid:90) ( x ,x ,x ) ∈ Ω ξ ( x , x , x ) dx dx dx = (cid:90) r = R out r = R in (cid:90) θ =1 θ =0 (cid:90) θ = N θ =0 (cid:101) ξ ( r, θ , θ ) 2 πτ r dr dθ dθ , (D.18)with (cid:101) ξ ( r, θ , θ ) = ξ ( x ( r, θ , θ ) , x ( r, θ , θ ) , x ( r, θ , θ )) .— Acknowledgement
Some of the theoretical and computational foundations of this work were laid out whileASB was a graduate student at the University of Minnesota, and later a postdoctoralfellow at the Lawrence Berkeley National Laboratory. ASB would like to acknowledgesupport of the Scientific Discovery through Advanced Computing (SciDAC) programfunded by U.S. Department of Energy, Office of Science, Advanced Scientific ComputingResearch and Basic Energy Sciences during his time at the Berkeley Lab, as well as the59upport of the following grants while at the University of Minnesota: AFOSR FA9550-15-1-0207, NSF-PIRE OISE-0967140, ONR N00014-14-1-0714 and the MURI projectFA9550-12-1-0458 (administered by AFOSR).ASB would like to acknowledge informative discussions and email communicationswith Phanish Suryanarayana (Georgia Institute of Technology), Carlos Garcia Cervera(Univ. of California, Santa Barbara), Eric Cances (Ecole des Ponts ParisTech), RichardJames (Univ. of Minnesota), Ryan Elliott (Univ. of Minnesota), Paul Garrett (Univ. ofMinnesota), Kaushik Bhattacharya (Caltech), Lin Lin (Univ. of California, Berkeley)and Chao Yang (Lawrence Berkeley National Lab). ASB would also like to thank NehaBairoliya (Univ. of Southern California) for her help in preparing some of the figures, andalso for providing encouragement and support during the preparation of this manuscript.Help from Hsuan Ming Yu (UCLA) in generating some of the simulation data in thepaper is also gratefully acknowledged.ASB would like to thank the Minnesota Supercomputing Institute (MSI) and UCLA’sInstitute for Digital Research and Education (IDRE) for making available the computingresources used in this work. —
References [1] B. Bhushan, Springer handbook of nanotechnology, Springer, 2017.[2] G. Cao, Nanostructures and nanomaterials: synthesis, properties and applications, World scientific,2004.[3] A. K. Geim, Random walk to graphene, , Nobel Lecture (2010).[4] R. E. Smalley, Dicsovering the fullerenes, , Nobel Lecture (1996).[5] M. Fulekar, Nanotechnology: Importance and Applications, 1st Edition, I K International Pub-lishing House, 2010.[6] Nanotechnology and you: Benefits and applications, , official website of the United States Na-tional Nanotechnology Initiative.[7] R. D. James, Objective structures, Journal of the Mechanics and Physics of Solids 54 (11) (2006)2354–2390.[8] A. S. Banerjee, Density functional methods for Objective Structures: Theory and simulationschemes, Ph.D. thesis, University of Minnesota, Minneapolis (2013).[9] T. Hahn, International tables for crystallography: volumes A,E, in: V. Kopsky, D. B. Litvin (Eds.),International Tables for Crystallography., Kluwer Academic Publishers, 2003.[10] K. Dayal, R. S. Elliott, R. D. James, Objective formulas, (In preparation) (2020).[11] Z. Ren, P.-X. Gao, A review of helical nanostructures: growth theories, synthesis strategies andproperties, Nanoscale 6 (16) (2014) 9366–9400.[12] David Goodsell, Molecule of the month: Amyloids (2015). doi:10.2210/rcsb_pdb/mom_2015_9 .[13] H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov,P. E. Bourne, The Protein Data Bank, Nucleic Acids Research 28 (1) (2000) 235–242. doi:10.1093/nar/28.1.235 .URL https://doi.org/10.1093/nar/28.1.235 [14] D. S. Goodsell, S. Dutta, C. Zardecki, M. Voigt, H. M. Berman, S. K. Burley, The rcsb pdbâĂIJmolecule of the monthâĂİ: inspiring a molecular view of biology, PLoS Biol 13 (5) (2015)e1002140.[15] T. Shimada, K. Nagano, J.-S. Wang, J. Wang, T. Kitamura, et al., Polar superhelices in ferroelectricchiral nanosprings, Scientific reports 6 (1) (2016) 1–9.[16] R. Naaman, D. H. Waldeck, Spintronics and chirality: Spin selectivity in electron transport throughchiral molecules, Annual review of physical chemistry 66 (2015) 263–281.
17] S. Dalum, P. Hedegård, Theory of chiral induced spin selectivity, Nano letters 19 (8) (2019) 5253–5259.[18] E. Medina, L. A. González-Arraga, D. Finkelstein-Shapiro, B. Berche, V. Mujica, Continuummodel for chiral induced spin selectivity in helical molecules, The Journal of Chemical Physics142 (19) (2015) 194308.[19] F. Giustino, Materials modelling using density functional theory: properties and predictions, Ox-ford University Press (UK), 2014.[20] J. Hafner, Ab-initio simulations of materials using vasp: Density-functional theory and beyond,Journal of computational chemistry 29 (13) (2008) 2044–2078.[21] T. Ziegler, Approximate density functional theory as a practical tool in molecular energetics anddynamics, Chemical Reviews 91 (5) (1991) 651–667.[22] A. Jain, Y. Shin, K. A. Persson, Computational predictions of energy materials using densityfunctional theory, Nature Reviews Materials 1 (1) (2016) 1–13.[23] J. Hafner, C. Wolverton, G. Ceder, Toward computational materials design: the impact of densityfunctional theory on materials research, MRS bulletin 31 (9) (2006) 659–668.[24] W. Kohn, L. J. Sham, Self-consistent equations including exchange and correlation effects, PhysicalReview 140 (4A) (1965) 1133–1138.[25] C. Le Bris (Ed.), Computational Chemistry, Vol. X of Handbook of Numerical Analysis, North-Holland, 2003.[26] G. Kresse, J. Furthmuller, Efficient iterative schemes for ab initio total-energy calculations usinga plane-wave basis set, Physical Review B 54 (1996) 11169–11186.[27] M. D. Segall, P. J. D. Lindan, M. J. Probert, C. J. Pickard, P. J. Hasnip, S. J. Clark, M. C.Payne, First-principles simulation: ideas, illustrations and the CASTEP code, Journal of Physics:Condensed Matter 14 (11) (2002) 2717.[28] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti,M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Ger-stmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Maz-zarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular andopen-source software project for quantum simulations of materials, Journal of Physics: CondensedMatter 21 (39) (2009).[29] X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Ver-straete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, P. Ghosez, J.-Y. Raty, D. Allan,First-principles computation of material properties: the ABINIT software project, ComputationalMaterials Science 25 (3) (2002) 478 – 492.[30] J. R. Chelikowsky, N. Troullier, K. Wu, Y. Saad, Higher order finite difference pseudopotentialmethod: An application to diatomic molecules, Phys. Rev. B 50 (1994) 11355–11364.[31] A. Castro, H. Appel, M. Oliveira, C. Rozzi, X. Andrade, F. Lorenzen, M. Marques, E. Gross,A. Rubio, Octopus: A tool for the application of time-dependent density functional theory, PhysicaStatus Solidi (B) Basic Research 243 (11) (2006) 2465–2488.[32] S. Ghosh, P. Suryanarayana, SPARC: Accurate and efficient finite-difference formulation and par-allel implementation of density functional theory: Isolated clusters, Computer Physics Communi-cations 212 (2017) 189–204.[33] S. Ghosh, P. Suryanarayana, SPARC: Accurate and efficient finite-difference formulation and par-allel implementation of density functional theory: Extended systems, Computer Physics Commu-nications 216 (2017) 109–125.[34] J. Pask, P. Sterne, Finite element methods in ab initio electronic structure calculations, Modellingand Simulation in Materials Science and Engineering 13 (3) (2005) R71.[35] J. Pask, B. Klein, P. Sterne, C. Fong, Finite-element methods in electronic-structure theory, Com-puter Physics Communications 135 (1) (2001) 1–34.[36] P. Motamarri, M. Nowak, K. Leiter, J. Knap, V. Gavini, Higher-order adaptive finite-elementmethods for Kohn-Sham density functional theory, Journal of Computational Physics 253 (2013)308–343.[37] N. W. Ashcroft, N. D. Mermin, Solid State Physics, 1st Edition, Brooks Cole, 1976.[38] C. Kittel, Introduction to Solid State Physics, eighth Edition, Wiley, 2004.[39] A. S. Banerjee, P. Suryanarayana, Cyclic density functional theory: A route to the first principles imulation of bending in nanostructures, Journal of the Mechanics and Physics of Solids 96 (2016)605–631.[40] A. O. Barut, R. Raczka, Theory of Group Representations and Applications, second revised Edi-tion, World Scientific Publishing Company, 1986.[41] G. B. Folland, A Course in Abstract Harmonic Analysis, 1st Edition, Studies in Advanced Math-ematics, Taylor & Francis, 1994.[42] C. Le Bris, P.-L. Lions, From atoms to crystals: a mathematical journey, Bulletin of the AmericanMathematical Society 42 (3) (2005) 291–363.[43] M. Defranceschi, C. Le Bris, Computing a molecule: A mathematical viewpoint, Journal of Math-ematical Chemistry 21 (1) (1997) 1–30.[44] A. S. Banerjee, R. S. Elliott, R. D. James, A spectral scheme for kohn–sham density functionaltheory of clusters, Journal of Computational Physics 287 (2015) 226–253.[45] S. Agarwal, A. S. Banerjee, On the use of helical waves for electronic structure calculations ofquasi-one-dimensional materials, (in preparation) (2020).[46] Y. Wei, B. Wang, J. Wu, R. Yang, M. L. Dunn, Bending rigidity and gaussian bending stiffness ofsingle-layered graphene, Nano letters 13 (1) (2012) 26–30.[47] I. Naumov, A. Bratkovsky, Gap opening in graphene by simple periodic inhomogeneous strain,Physical Review B 84 (24) (2011) 245444.[48] T. Dumitrica, R. D. James, Objective molecular dynamics, Journal of the Mechanics and Physicsof Solids 55 (10) (2007) 2206 – 2236.[49] D. Zhang, E. Akatyeva, T. Dumitrică, Bending ultrathin graphene at the margins of continuummechanics, Phys. Rev. Lett. 106 (2011) 255503.[50] T. Dumitrica, Computational nanomechanics of quasi-one-dimensional structures in a symmetry-adapted tight binding framework, in: V. Barsan, A. Aldea (Eds.), Trends in Nanophysics, Vol. 0of Engineering Materials, Springer Berlin Heidelberg, 2010, pp. 29–55.[51] J. Ma, Y. Ni, S. Volz, T. Dumitrică, Thermal transport in single-walled carbon nanotubes underpure bending, Physical Review Applied 3 (2) (2015) 024014.[52] P. Koskinen, O. O. Kit, Efficient approach for simulating distorted materials, Physical reviewletters 105 (10) (2010) 106401.[53] P. Koskinen, Electronic and optical properties of carbon nanotubes under pure bending, Phys.Rev. B 82 (2010) 193409.[54] P. Koskinen, Graphene nanoribbons subject to gentle bends, Physical Review B 85 (20) (2012)205429.[55] O. O. Kit, L. Pastewka, P. Koskinen, Revised periodic boundary conditions: Fundamentals, elec-trostatics, and the tight-binding approximation, Physical Review B 84 (15) (2011) 155431.[56] S. Ismail-Beigi, T. Arias, Ab initio study of screw dislocations in Mo and Ta: a new picture ofplasticity in bcc transition metals, Physical Review Letters 84 (7) (2000) 1499.[57] J. A. Hauch, D. Holland, M. Marder, H. L. Swinney, Dynamic fracture in single crystal silicon,Physical Review Letters 82 (19) (1999) 3823.[58] G. Cocco, E. Cadelano, L. Colombo, Gap opening in graphene by shear strain, Physical Review B81 (24) (2010) 241412.[59] P. Koskinen, V. Mäkinen, Density-functional tight-binding for beginners, Computational MaterialsScience 47 (1) (2009) 237–253.[60] S. Ghosh, A. S. Banerjee, P. Suryanarayana, Symmetry-adapted real-space density functionaltheory for cylindrical geometries: Application to large group-IV nanotubes, Physical Review B100 (12) (2019) 125143.[61] P. DâĂŹArco, Y. Noel, R. Demichelis, R. Dovesi, Single-layered chrysotile nanotubes: A quantummechanical ab initio simulation, The Journal of chemical physics 131 (20) (2009) 204701.[62] R. Dovesi, V. Saunders, C. Roetti, R. Orlando, C. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll,N. Harrison, I. Bush, et al., Crystal17 (2017).[63] C. T. White, D. H. Robertson, J. W. Mintmire, Helical and rotational symmetries of nanoscalegraphitic tubules, Physical Review B 47 (1993) 5485–5488.[64] E. Chang, G. Bussi, A. Ruini, E. Molinari, First-principles approach for the calculation of opticalproperties of one-dimensional systems with helical symmetry: The case of carbon nanotubes,Physical Review B 72 (2005) 195423.[65] Y. Noel, P. D’arco, R. Demichelis, C. M. Zicovich-Wilson, R. Dovesi, On the use of symmetry n the ab initio quantum mechanical simulation of nanotubes and related materials, Journal ofComputational Chemistry 31 (4) (2010) 855–862.[66] E. Chang, G. Bussi, A. Ruini, E. Molinari, Excitons in carbon nanotubes: An Ab Initio symmetry-based approach, Physical Review Letters 92 (2004) 196401.[67] F. Feng, P. Plucinsky, R. D. James, Phase transformations and compatibility in helical structures,Journal of the Mechanics and Physics of Solids (2019).[68] F. Feng, Phase transformation in helical structures: Theory and application, Ph.D. thesis, Univer-sity of Minnesota (2018).[69] W. Miller, Symmetry Groups and Their Applications, Vol. 50 of Pure and Applied MathematicsSeries, Elsevier Science & Technology Books, 1972.[70] H. M. Yu, A. S. Banerjee, Symmetry-adapted real-space density functional theory in helical coor-dinates: Applications to torsional deformations in group-IV nanotubes, (in preparation) (2020).[71] I. Nikiforov, B. Hourahine, B. Aradi, T. Frauenheim, T. Dumitrică, Ewald summation on a helix: Aroute to self-consistent charge density-functional based tight-binding objective molecular dynamics,The Journal of Chemical Physics 139 (9) (2013) 094110.[72] M. Defranceschi, C. Le Bris (Eds.), Mathematical Models and Methods for Ab Initio QuantumChemistry, Vol. 74 of Lecture Notes in Chemistry, Springer, 2000.[73] R. M. Martin, Electronic Structure: Basic Theory and Practical Methods, 1st Edition, CambridgeUniversity Press, 2004.[74] F. Bloch, Über die quantenmechanik der elektronen in kristallgittern, Zeitschrift für Physik AHadrons and Nuclei 52 (1929) 555–600.[75] P. Koskinen, O. O. Kit, Efficient approach for simulating distorted materials, Phys. Rev. Lett. 105(2010) 106401.[76] O. O. Kit, L. Pastewka, P. Koskinen, Revised periodic boundary conditions: Fundamentals, elec-trostatics, and the tight-binding approximation, Phys. Rev. B 84 (2011) 155431.[77] I. Milosević, B. Dakić, M. Damnjanović, Generalized bloch states and potentials of nanotubes andother quasi-1d systems, Journal of Physics A: Mathematical and General 39 (38) (2006) 11833.[78] B. Dakić, M. Damnjanović, I. Milosević, Generalized bloch states and potentials of nanotubesand other quasi-1d systems ii, Journal of Physics A: Mathematical and Theoretical 42 (12) (2009)125202.[79] J. Mintmire, C. White, First-principles band structures of armchair nanotubes, Applied PhysicsA: Materials Science and Processing 67 (1998) 65–69.[80] J. Mintmire, C. White, Electronic and structural properties of carbon nanotubes, Carbon 33 (7)(1995) 893 – 902.[81] E. B. Barros, A. Jorio, G. G. Samsonidze, R. B. Capaz, A. G. S. Filho, J. M. Filho, G. Dressel-haus, M. S. Dresselhaus, Review on the symmetry-related properties of carbon nanotubes, PhysicsReports 431 (6) (2006) 261 – 302.[82] F. Odeh, J. B. Keller, Partial differential equations with periodic coefficients and bloch waves incrystals, Journal of Mathematical Physics 5 (1964) 1499–1504.[83] W. Dörfler, A. Lechleiter, M. Plum, G. Schneider, C. Wieners, Photonic crystals: Mathematicalanalysis and numerical approximation, Vol. 42, Springer Science & Business Media, 2011.[84] C. Wilcox, Theory of Bloch waves, Journal d’Analyse MathÃľmatique 33 (1978) 146–167.[85] M. Reed, B. Simon, Analysis of Operators, Vol. IV of Methods of Modern Mathematical Physics,Academic Press, 1978.[86] M. Hoffmann-Ostenhof, T. Hoffmann-Ostenhof, R. Ahlrichs, J. Morgan, On the exponential falloff of wavefunctions and electron densities, in: K. Osterwalder (Ed.), Mathematical Problems inTheoretical Physics, Vol. 116 of Lecture Notes in Physics, Springer Berlin / Heidelberg, 1980, pp.62–67.[87] R. Ahlrichs, M. Hoffmann-Ostenhof, T. Hoffmann-Ostenhof, J. D. Morgan, Bounds on the decayof electron densities with screening, Physical Review A 23 (5) (1981) 2106–2117.[88] J. Han, M. L. Tiago, T.-L. Chan, J. R. Chelikowsky, Real space method for the electronic structureof one-dimensional periodic systems, The Journal of chemical physics 129 (14) (2008) 144109.[89] L. C. Evans, Partial Differential Equations, Vol. 19 of Graduate Studies in Mathematics, AmericanMathematical Society, 1998.[90] M. Renardy, R. C. Rogers, An Introduction to Partial Differential Equations, 2nd Edition, Vol. 13of Texts in Applied Mathematics, Springer, 2004.
91] T. Kato, Perturbation Theory for Linear Operators, Classics in Mathematics, Springer, 1995.[92] F. Rellich, J. Berkowitz, Perturbation theory of eigenvalue problems, CRC Press, 1969.[93] R. McWeeny, Symmetry:An Introduction to Group Theory and Its Applications, 1st Edition,Dover, 2002.[94] G. B. Folland, Real Analysis: Modern Techniques and Their Applications, 2nd Edition, Wiley,1999.[95] D. Jüstel, The zak transform on strongly proper g-spaces and its applications, Journal of theLondon Mathematical Society 97 (1) (2018) 47–76.[96] D. Jüstel, Radiation for the analysis of molecular structures with non-crystalline symmetry: mod-elling and representation theoretic design, Ph.D. thesis, Technische Universität München (2014).[97] P. D. Hislop, I. M. Sigal, Introduction to spectral theory: With applications to Schrödinger oper-ators, Vol. 113, Springer Science & Business Media, 2012.[98] I. Catto, C. Le Bris, P.-L. Lions, On the thermodynamic limit for hartree-fock type models, in:Annales de l’IHP Analyse non linéaire, Vol. 18, 2001, pp. 687–760.[99] E. Cances, V. Ehrlacher, D. Gontier, A. Levitt, D. Lombardi, Numerical quadrature in the brillouinzone for periodic schrödinger operators, Numerische Mathematik 144 (3) (2020) 479–526.[100] E. Cancès, A. Deleurence, M. Lewin, A new approach to the modeling of local defects in crystals:the reduced hartree-fock case, Communications in Mathematical Physics 281 (1) (2008) 129–177.[101] S. Goedecker, Linear scaling electronic structure methods, Reviews of Modern Physics 71 (4) (1999)1085.[102] A. Anantharaman, E. Cancès, Existence of minimizers for kohn–sham models in quantum chem-istry, in: Annales de l’Institut Henri Poincare (C) Non Linear Analysis, Vol. 26, Elsevier, 2009,pp. 2425–2455.[103] X.-C. Wang, T. Blesgen, K. Bhattacharya, M. Ortiz, A variational framework for spectral approx-imations of kohn–sham density functional theory, Archive for Rational Mechanics and Analysis221 (2) (2016) 1035–1075.[104] P. Garrett, A good spectral theorem, .[105] K. Schmüdgen, Unbounded self-adjoint operators on Hilbert space, Vol. 265, Springer Science &Business Media, 2012.[106] G. Teschl, Mathematical Methods in Quantum Mechanics, Graduate Studies in Mathematics,American Mathematical Society, 2009.[107] M. Benzi, P. Boito, N. Razouk, Decay properties of spectral projectors with applications to elec-tronic structure, SIAM review 55 (1) (2013) 3–64.[108] N. Troullier, J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Physical reviewB 43 (3) (1991) 1993.[109] J. R. Chelikowsky, Introductory Quantum Mechanics with MATLAB: For Atoms, Molecules, Clus-ters, and Nanocrystals, John Wiley & Sons, 2019.[110] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys.Rev. Lett. 77 (1996) 3865–3868.[111] L. Kleinman, D. Bylander, Efficacious form for model pseudopotentials, Physical Review Letters48 (20) (1982) 1425.[112] N. Troullier, J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Physical ReviewB 43 (3) (1991) 1993.[113] J. E. Pask, P. A. Sterne, Real-space formulation of the electrostatic potential and total energy ofsolids, Phys. Rev. B 71 (2005) 113101.[114] P. Suryanarayana, D. Phanish, Augmented lagrangian formulation of orbital-free density functionaltheory, Journal of Computational Physics 275 (2014) 524–538.[115] P. Motamarri, M. Iyer, J. Knap, V. Gavini, Higher-order adaptive finite-element methods fororbital-free density functional theory, Journal of Computational Physics 231 (20) (2012) 6596–6621.[116] P. Motamarri, V. Gavini, Configurational forces in electronic structure calculations using kohn-sham density functional theory, Physical Review B 97 (16) (2018) 165132.[117] J. Harris, Simplified method for calculating the energy of weakly interacting fragments, PhysicalReview B 31 (4) (1985) 1770.[118] W. M. C. Foulkes, R. Haydock, Tight-binding models and density-functional theory, Physical eview B 39 (17) (1989) 12520.[119] H. Kikuji, O. Tomoya, F. Yoshitaka, T. Shigeru, First-principles calculations in real-space formal-ism: electronic configurations and transport properties of nanostructures, World Scientific, 2005.[120] H. J. Monkhorst, J. D. Pack, Special points for brillouin-zone integrations, Physical review B13 (12) (1976) 5188.[121] I. I. Geru, Time-Reversal Symmetry, Springer, 2018.[122] A. S. Banerjee, L. Lin, P. Suryanarayana, C. Yang, J. E. Pask, Two-level chebyshev filter basedcomplementary subspace method: pushing the envelope of large-scale electronic structure calcula-tions, Journal of chemical theory and computation 14 (6) (2018) 2930–2946.[123] D. A. Mazziotti, Spectral difference methods for solving differential equations, Chemical physicsletters 299 (5) (1999) 473–480.[124] Y. Zhou, J. R. Chelikowsky, Y. Saad, Chebyshev-filtered subspace iteration method free of sparsediagonalization for solving the kohn–sham equation, Journal of Computational Physics 274 (2014)770–782.[125] Y. Zhou, Y. Saad, M. L. Tiago, J. R. Chelikowsky, Self-consistent-field calculations usingChebyshev-filtered subspace iteration, Journal of Computational Physics 219 (2006) 172–184.[126] Y. Zhou, Y. Saad, M. L. Tiago, J. R. Chelikowsky, Parallel self-consistent-field calculations viaChebyshev-filtered subspace acceleration, Phys. Rev. E 74 (2006) 066704.[127] F. Gygi, G. Galli, Real-space adaptive-coordinate electronic-structure calculations, Physical Re-view B 52 (4) (1995) R2229.[128] Y. Saad, M. H. Schultz, Gmres: A generalized minimal residual algorithm for solving nonsymmetriclinear systems, SIAM Journal on scientific and statistical computing 7 (3) (1986) 856–869.[129] Y. Saad, Iterative methods for sparse linear systems, 2nd Edition, SIAM, 2003.[130] MATLAB, 9.7.0.1190202 (R2019b), the MathWorks Inc., Natick, Massachusetts (2019).[131] A. S. Banerjee, P. Suryanarayana, J. E. Pask, Periodic Pulay method for robust and efficientconvergence acceleration of self-consistent field iterations, Chemical Physics Letters 647 (2016)31–35.[132] P. Pulay, Convergence acceleration of iterative sequences. the case of scf iteration, Chemical PhysicsLetters 73 (2) (1980) 393–398.[133] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumbsch, Structural relaxation made simple,Physical Review Letters 97 (17) (2006) 170201.[134] J. P. Perdew, Y. Wang, Accurate and simple analytic representation of the electron-gas correlationenergy, Physical Review B 45 (1992) 13244–13249.[135] R. A. Èvarestov, Theoretical modeling of inorganic nanostructures: Symmetry and ab-initio cal-culations of nanolayers, nanotubes and nanowires, Springer, 2015.[136] E. S. Reich, et al., Phosphorene excites materials scientists, Nature 506 (7486) (2014) 19.[137] S. Das, W. Zhang, M. Demarteau, A. Hoffmann, M. Dubey, A. Roelofs, Tunable transport gap inphosphorene, Nano letters 14 (10) (2014) 5733–5739.[138] A. Carvalho, M. Wang, X. Zhu, A. S. Rodin, H. Su, A. H. C. Neto, Phosphorene: from theory toapplications, Nature Reviews Materials 1 (11) (2016) 1–16.[139] H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, P. D. Ye, Phosphorene: an unexplored2d semiconductor with a high hole mobility, ACS nano 8 (4) (2014) 4033–4041.[140] L. Kou, C. Chen, S. C. Smith, Phosphorene: fabrication, properties, and applications, The journalof physical chemistry letters 6 (14) (2015) 2794–2805.[141] A. Rodin, A. Carvalho, A. C. Neto, Strain-induced gap modification in black phosphorus, PhysicalReview Letters 112 (17) (2014) 176801.[142] J. Guan, Z. Zhu, D. Tománek, Phase coexistence and metal-insulator transition in few-layer phos-phorene: a computational study, Physical Review Letters 113 (4) (2014) 046804.[143] H. Guo, N. Lu, J. Dai, X. Wu, X. C. Zeng, Phosphorene nanoribbons, phosphorus nanotubes, andvan der waals multilayers, The Journal of Physical Chemistry C 118 (25) (2014) 14051–14059.[144] V.-T. Nguyen, M.-Q. Le, Atomistic simulation of the uniaxial compression of black phosphorenenanotubes, Vietnam Journal of Mechanics 40 (3) (2018) 243–50.[145] J. Cao, J. Shi, Y. Hu, M. Wu, C. Ouyang, B. Xu, Lithium ion adsorption and diffusion on blackphosphorene nanotube: A first-principles study, Applied Surface Science 392 (2017) 88–94.[146] V. Sorkin, Y. Zhang, Mechanical properties of phosphorene nanotubes: a density functional tight-binding study, Nanotechnology 27 (39) (2016) 395701.199] G. De Nittis, M. Lein, Exponentially localized wannier functions in periodic zero flux magneticfields, Journal of mathematical physics 52 (11) (2011) 112103.[200] G. Panati, C. Sparber, S. Teufel, Geometric currents in piezoelectricity, Archive for rational me-chanics and analysis 191 (3) (2009) 387–422.