Symmetry-adapted real-space density functional theory for cylindrical geometries: application to large X (X=C, Si, Ge, Sn) nanotubes

TL;DR

This paper introduces a symmetry-adapted real-space density functional theory approach for cylindrical geometries, enabling efficient large-scale simulations of nanotubes and Xene sheets with detailed analysis of their electronic and mechanical properties.

Contribution

It develops a high-order finite-difference implementation incorporating symmetry reduction, allowing accurate and efficient ab initio calculations of large nanotubes and related structures.

Findings

Nanotubes are semiconducting with bandgap inversely proportional to radius.

Bending moduli of Xene sheets show Kirchhoff-Love behavior and anisotropy.

The method scales well and enables simulations of micron-sized nanotubes.

Abstract

We present a symmetry-adapted real-space formulation of Kohn-Sham density functional theory for cylindrical geometries and apply it to the study of large X (X=C, Si, Ge, Sn) nanotubes. Specifically, starting from the Kohn-Sham equations posed on all of space, we reduce the problem to the fundamental domain by incorporating cyclic and periodic symmetries present in the angular and axial directions of the cylinder, respectively. We develop a high-order finite-difference parallel implementation of this formulation, and verify its accuracy against established planewave and real-space codes. Using this implementation, we study the band structure and bending properties of X nanotubes and Xene sheets, respectively. Specifically, we first show that zigzag and armchair X nanotubes with radii in the range 1 to 5 nm are semiconducting. In particular, we find an inverse linear dependence of the bandgap with respect to the radius for all nanotubes, other than the armchair and zigzag type III carbon variants, for which we find an inverse quadratic dependence. Next, we exploit the connection between cyclic symmetry and uniform bending deformations to calculate the bending moduli of Xene sheets in both zigzag and armchair directions. We find Kirchhoff-Love type bending behavior for all sheets, with graphene and stanene possessing the largest and smallest moduli, respectively. In addition, other than graphene, the sheets demonstrate significant anisotropy, with larger bending moduli along the armchair direction. Finally, we demonstrate that the proposed approach has very good parallel scaling and is highly efficient, enabling ab initio simulations of unprecedented size for systems with a high degree of cyclic symmetry. In particular, we show that even micron-sized nanotubes can be simulated with modest computational effort.

Tables

Table 1. Table 1: Equilibrium lattice parameters for the Xene sheets. The numbers in parenthesis are values from literature.

Xene	$a$ (Å)	$\delta$ (Å)
C	1.407 (1.408Kerszberg and Suryanarayana (2015))	-
Si	2.200 (2.207Buda et al. (2017))	0.404 (0.437Buda et al. (2017))
Ge	2.232 (2.290Buda et al. (2017))	0.566 (0.647Buda et al. (2017))
Sn	2.522 (2.611 Zhou et al. (2016))	0.699 (0.822 Zhou et al. (2016))

Table 2. Table 2: Real- and η 𝜂 \eta -space discretization parameters for the X nanotube simulations in Sections  V.2 and V.3 , from which the physical properties of interest, i.e., variation of bandgap with radius of X nanotubes and the bending moduli of Xene sheets, are calculated.

X	$h$ (Bohr)	$N_{\eta }$
C	0.125	Armchair: $21$, Zigzag: $13$
Si	0.250	Armchair: $15$, Zigzag: $9$
Ge	0.200	Armchair: $15$, Zigzag: $9$
Sn	0.250	Armchair: $15$, Zigzag: $9$

Table 3. Table 3: Power law exponents for the decay of bandgap with radius in X nanotubes.

X	Armchair		Zigzag
		Type I	Type II	Type III
C	-2.01	-0.96	-1.03	-1.95
Si	-1.00	-1.00	-1.03	-1.03
Ge	-1.01	-1.03	-1.04	-1.04
Sn	-1.02	-1.05	-1.06	-1.14

Table 4. Table 4: Bending modulus ( D 𝐷 D ) for the Xene sheets in the armchair and zigzag directions.

Xene	$D$ (eV)
	Armchair	Zigzag
C	1.57	1.50
Si	0.41	0.54
Ge	0.38	0.56
Sn	0.26	0.46

Equations

$$\bar{\mathcal{F}}(\Psi ,\mathbf{g},\text{X})={\bar{T}}_s(\Psi ,\mathbf{g})$$

$$\rho (\mathbf{r})=2\sum _{n=1}^{{\bar{N}}_s}{g}_n|{\psi }_n(\mathbf{r}){|}^2,$$

$${\bar{T}}_s(\Psi ,\mathbf{g})=-\sum _{n=1}^{{\bar{N}}_s}{g}_n{\int }_{{\mathbb{R}}^3}{\psi }_n^{\ast }(\mathbf{r}){\nabla }^2{\psi }_n(\mathbf{r})\mathrm{d}\mathbf{r}.$$

$$E_{xc}(\rho )={\int }_{{\mathbb{R}}^3}{\epsilon }_{xc}(\rho (\mathbf{r}))\rho (\mathbf{r})\mathrm{d}\mathbf{r}.$$

$$\bar{K}(\Psi ,\mathbf{g},\text{X})=2\sum _{n=1}^{{\bar{N}}_s}{g}_n\sum _J\sum _{p\in {\mathcal{A}}_J}{\gamma }_{J;p}{\left|{\int }_{{\mathbb{R}}^3}{\chi }_{J;p}^{\ast }({\text{x}}_J,\text{r}){\psi }_n(\text{r})\mathrm{d}\text{r}\right|}^2.$$

$$\displaystyle\bar{E}_{el}(\rho,\textbf{X})=\max_{\phi}\bigg{\{}-\frac{1}{8\pi}\int_{\mathbb{R}^{3}}|\nabla\phi({\mathbf{r}})|^{2}\,\mathrm{d{\mathbf{r}}}+\int_{\mathbb{R}^{3}}\big{(}\rho({\mathbf{r}})+b({\mathbf{r}},\textbf{X})\big{)}\phi({\mathbf{r}})\,\mathrm{d{\mathbf{r}}}\bigg{\}}+\bar{E}_{sc}(\textbf{X})\,.$$

$$\displaystyle\bar{S}({\mathbf{g}})=-2k_{B}\sum_{n=1}^{\bar{N}_{s}}\big{(}g_{n}\log g_{n}+(1-g_{n})\log(1-g_{n})\big{)}\,,$$

$${\bar{\mathcal{F}}}_0(\text{X})$$

$$s.t.{\int }_{{\mathbb{R}}^3}{\psi }_i(\mathbf{x}){\psi }_j(\mathbf{x})\mathrm{d}\mathbf{x}$$

$$\left(\bar{\mathcal{H}}(\Psi ,\mathbf{g},\text{X})\equiv -\frac{1}{2}{\nabla }^2+V_{xc}+\varphi +{\bar{\mathcal{V}}}_{nl}\right){\psi }_n$$

$$-\frac{1}{4\pi }{\nabla }^2\varphi (\mathbf{r},\text{X})=\rho (\mathbf{r})+b(\mathbf{r},\text{X}),$$

$$\displaystyle g_{n}=\bigg{(}1+\exp{\bigg{(}\frac{\lambda_{n}-\lambda_{F}}{k_{B}T}\bigg{)}}\bigg{)}^{-1}\,,$$

$$[{\bar{\mathcal{V}}}_{nl}f](\mathbf{r})=\sum _J\sum _{p\in {\mathcal{A}}_J}{\gamma }_{J;p}{\chi }_{J;p}({\mathbf{x}}_J,\mathbf{r}){\int }_{{\mathbb{R}}^3}{\chi }_{\mathbf{x};p}^{\ast }(\mathbf{x},\text{y})f(\text{y})\mathrm{d}\text{y}.$$

$$\displaystyle\bar{\mathcal{F}}_{0}(\textbf{X})=2\sum_{n=1}^{\bar{N}_{s}}g_{n}\lambda_{n}+\bar{E}_{xc}(\rho)-\int_{{\mathbb{R}}^{3}}V_{xc}(\rho({\mathbf{r}}))\rho({\mathbf{r}})\,\mathrm{d{\mathbf{r}}}+\frac{1}{2}\int_{{\mathbb{R}}^{3}}\big{(}b({\mathbf{r}},\mathbf{X})-\rho({\mathbf{r}})\big{)}\phi({\mathbf{r}},\mathbf{X})\,\mathrm{d{\mathbf{r}}}+\bar{E}_{sc}(\mathbf{X})-T\bar{S}({\mathbf{g}})\,,$$

$${\mathbf{f}}_J$$

$$\Omega:=\left\{(r,\theta,z)\in{\mathbb{R}}^{3}\big{|}R_{in}\leq r\leq R_{out},0\leq\theta\leq\Theta\leq 2\pi,0\leq z\leq H\right\}\,,$$

$$\partial \Omega =\partial R_{in}\bigcup \partial R_{out}\bigcup \partial {\vartheta }_0\bigcup \partial {\vartheta }_{\Theta }\bigcup \partial {\mathcal{Z}}_0\bigcup \partial {\mathcal{Z}}_H,$$

$$\mathfrak{T}=\{{\text{t}}_{\mu }:\mu \in \mathbb{Z}\},\text{with}{\text{t}}_{\mu }=\mu H{\text{e}}_z.$$

$$\mathfrak{C}=\{{\mathfrak{R}}^{\zeta }:\zeta =0,1,2,\dots ,\mathfrak{N}-1\},\text{with}{\mathfrak{R}}^{\zeta }=\left(\begin{array}{ccc}\cos \frac{2\pi \zeta }{\mathfrak{N}}& -\sin \frac{2\pi \zeta }{\mathfrak{N}}& 0\text{vspace}0.05in\\ \sin \frac{2\pi \zeta }{\mathfrak{N}}& \cos \frac{2\pi \zeta }{\mathfrak{N}}& 0\text{vspace}0.05in\\ 0& 0& 1\end{array}\right).$$

$$\mathcal{G}=\{{{\rm Y}}_{\zeta ,\mu }=({\mathfrak{R}}^{\zeta }|{\text{t}}_{\mu }):{\mathfrak{R}}^{\zeta }\in \mathfrak{C},{\text{t}}_{\mu }\in \mathfrak{T}\},$$

$$\text{X}=\{{{\rm Y}}_{\zeta ,\mu }\circ {\mathbf{R}}_J={\mathfrak{R}}^{\zeta }{\mathbf{R}}_J+{\text{t}}_{\mu }\},\text{with}J=1,2,\dots ,N,\zeta =0,1,2,\dots ,\mathfrak{N}-1,\mu \in \mathbb{Z}.$$

$$\displaystyle\mathcal{G}\circ{\mathbf{R}}_{J}=\big{\{}\Upsilon_{\zeta,\mu}\circ{\mathbf{R}}_{J}:\Upsilon_{\zeta,\mu}\in\mathcal{G}\big{\}}\,.$$

$$\text{X}=\bigcup _{J=1}^N\mathcal{G}\circ {\mathbf{R}}_J.$$

$$\displaystyle\widehat{\mathcal{G}}=\bigg{\{}e^{2\pi i\big{(}\frac{\zeta}{\mathfrak{N}}{\nu}+\frac{\mu H}{2\pi}{\eta}\big{)}}:\,$$

$$\displaystyle\rho({\mathbf{r}})=2\sum_{n=1}^{N_{s}}\bigg{(}\frac{1}{\mathfrak{N}}\sum_{\nu=0}^{\mathfrak{N}-1}\fint g_{n}(\nu,\eta)\,|\psi_{n}({\mathbf{r}},\nu,\eta)|^{2}\,\mathrm{d\eta}\,\bigg{)}.$$

$$\displaystyle\psi_{n}(\Upsilon_{\zeta,\mu}^{-1}\circ{\mathbf{r}},\nu,\eta)=e^{2\pi i\big{(}\frac{\zeta}{\mathfrak{N}}{\nu}+\frac{\mu H}{2\pi}{\eta}\big{)}}\psi_{n}({\mathbf{r}},\nu,\eta)\,,$$

$$\displaystyle\psi_{n}(\Upsilon_{\zeta,\mu}\circ{\mathbf{r}},\nu,\eta)=e^{-2\pi i\big{(}\frac{\zeta}{\mathfrak{N}}{\nu}+\frac{\mu H}{2\pi}{\eta}\big{)}}\psi_{n}({\mathbf{r}},\nu,\eta)\,.$$

$${\psi }_n(r,\theta =\Theta ,z,\nu ,\eta )=e^{-\frac{2\pi i\nu }{\mathfrak{N}}}{\psi }_n(r,\theta =0,z,\nu ,\eta ),$$

$${\psi }_n(r,\theta ,z=H,\nu ,\eta )=e^{-i\eta H}{\psi }_n(r,\theta ,z=0,\nu ,\eta ).$$

$${\psi }_n(r=R_{in},\theta ,z,\nu ,\eta )={\psi }_n(r=R_{out},\theta ,z,\nu ,\eta )=0.$$

$$\varphi (r,\theta =\Theta ,z,\mathbf{R},\mathcal{G})=\varphi (r,\theta =0,z,\mathbf{R},\mathcal{G}),$$

Full text

I Introduction

Over the course of the past few decades, ab-initio calculations based on Kohn-Sham density functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) have become a mainstay of computational materials research due to their predictive power and ability to provide fundamental insights into materials properties and behavior. Indeed, a relatively large fraction of the computational resources worldwide are now devoted to first principles DFT calculations. The widespread popularity of DFT can be attributed to its generality, relative simplicity, and high accuracy-to-cost ratio when compared to other such ab-initio methods Burke (2012); Becke (2014). However, even though the cost of DFT is significantly less than the more accurate wavefunction-based theories, the efficient solution of the Kohn-Sham equations remains a challenging task, which severely restricts the range of physical systems that can be investigated. In particular, the computational cost and memory storage scale cubically and quadratically with respect to system size, respectively, and are typically associated with large prefactors. Moreover, the scalability in the context of high performance parallel computing suffers from the global nature of the orthonormality constraint posed on the Kohn-Sham orbitals.

The planewave pseudopotential method Martin (2004) is one of the most popular techniques for performing DFT calculations. In this approach, the Kohn-Sham equations are discretized using the Fourier basis, a complete and systematically improvable set in which convergence can be controlled by a single parameter. The planewave method is therefore not only accurate and simple to use, but is highly efficient on small to moderate computational resources due to the use of optimized Fast Fourier Transforms (FFTs). However, the planewave method suffers from a few limitations, including the following. The Fourier basis enforces periodic boundary conditions, whereby finite systems such as molecules and semi-infinite systems such as nanotubes require the introduction of artificial periodicity, with possibly large regions of vacuum between periodic replicas. This limitation also requires the introduction of an artificial neutralizing background density when treating charged systems. Furthermore, the global nature of the Fourier basis prevents the development of linear-scaling methods Goedecker (1999); Bowler and Miyazaki (2012). Finally, the reliance on FFTs hampers parallel scaling, thus limiting the length and time scales that can be reached.

In view of the aforementioned limitations of the planewave method, a number of alternate representations have been developed over the last two decades that are not only systematically improvable but also localized Arias (1999); Beck (2000); Pask and Sterne (2005a); Saad et al. (2010); Skylaris et al. (2005); Hernández et al. (1997); Suryanarayana et al. (2011); Lin et al. (2012). Among these approaches, real-space finite-difference methods Beck (2000); Saad et al. (2010)—for which computational locality is maximized by discretizing all quantities of interest on a real-space grid using high-order finite-differences—are some of the most mature and widely used to date. In these methods, convergence is again controlled by a single parameter, i.e., the mesh-size or grid spacing. In addition, any of Dirichlet, periodic, and Bloch-periodic boundary conditions (and combinations thereof) can be accommodated, whereby finite, semi-infinite, and charged systems, as well as bulk 3D systems can all be accurately and efficiently studied. Moreover, the locality of the discretization and decay of the density matrix within this representation Suryanarayana (2017) allows for the development of linear-scaling methods Osei-Kuffuor and Fattebert (2014); Suryanarayana (2013); Suryanarayana et al. (2017). Finally, large-scale parallel computational resources can be efficiently leveraged by virtue of the method’s simplicity, locality, and freedom from communication-intensive transforms such as FFTs. The limitations of real-space methods include the larger number of degrees of freedom/atom and the lack of effective preconditioners, when compared to the planewave method.

In recent years, there has been a significant increase in the efficiency of real-space methods due to a number of advances. Since early work in this area Bernholc et al. (1991); Chelikowsky et al. (1994); Briggs et al. (1995); Seitsonen et al. (1995), the degrees of freedom/atom required to obtain accurate ground state properties has been notably reduced by double-grid techniques Ono and Hirose (1999), ultrasoft pseudopotential formulations Hodak et al. (2007), projector augmented wave methods Mortensen et al. (2005), high-order integration Bobbitt et al. (2015), reformulation of nonlocal pseudopotential components Hirose et al. (2005); Sharma and Suryanarayana (2018a), and reduction of the eigenproblem by discontinuous projection Xu et al. (2018). Moreover, the need for effective preconditioners has been circumvented by substituting traditional iterative eigensolvers with the Chebyshev-polynomial filtered subspace iteration (CheFSI) Zhou et al. (2006a). These and other advances have made large systems containing thousands of atoms amenable to real-space finite-difference methods Alemany et al. (2008). In fact, they are now able to outperform established planewave codes in the context of both finite Ghosh and Suryanarayana (2017a) and extended Ghosh and Suryanarayana (2017b) systems. However, just like planewave methods, real-space methods have been restricted to affine (primarily Cartesian) coordinate systems. Though these are ideally suited for the various crystal systems, curvilinear coordinate systems provide an opportunity for the better description of systems with curved geometries.

1D nanostructures possessing a cylindrical-type geometry, e.g., nanotubes, nanowires, and nanorods, have received a lot of attention in the past three decades due to their unusual and fascinating material properties Martel et al. (1998); Javey et al. (2003); Popov (2004); Gong et al. (2009); Park et al. (2009); Wu et al. (2012); Park et al. (2011); Li et al. (2011); Zhao et al. (2006); Patzke et al. (2002); Law et al. (2004). This is also true for their 2D counterparts Xu et al. (2013); Bhimanapati et al. (2015); Butler et al. (2013); Naguib et al. (2014); Fiori et al. (2014); Koppens et al. (2014), which have risen to prominence after the discovery of graphene (Novoselov et al., 2004). Though these 2D structures are originally planar, they take up cylindrical-type geometries when subject to bending deformations, a common and technologically relevant mode of deformation in such materials Schniepp et al. (2008); Wei et al. (2012). These cylindrical-type geometries are indeed best described while working in a cylindrical coordinate system. In particular, such a choice enables the implementation to be compatible with the (possibly large) rotational/cyclic symmetry that could be present in the system. Indeed, such symmetry is commonly found in 1D nanostructures, and due to its connection with uniform bending deformations (James, 2006; Banerjee and Suryanarayana, 2016), also while studying 2D nanostructures subject to bending deformations. However, being restricted to affine coordinate systems, current DFT methods are unable to fully exploit the cyclic symmetry to reduce the computational cost, which is critical for 1D nanostructures such as nanotubes with large radii Shin et al. (2004); McGary et al. (2006); Macak et al. (2008) as well as for bending of 2D nanostructures with radii of curvature representative of those realized in experiments. The suitability of the real-space method for cylindrical coordinates and its attractive features outlined above provide the motivation for the current effort.

In this work, we present a symmetry-adapted real-space formulation of Kohn-Sham DFT for cylindrical geometries and apply it to the study of large X (X=C, Si, Ge, Sn) nanotubes. Specifically, we incorporate the cyclic and periodic symmetries present in the angular and axial directions of the cylinder, respectively, to reduce the Kohn-Sham equations that are originally posed on all of space to the fundamental domain. We develop a parallel implementation of this formulation using high-order finite-differences, and verify its accuracy by benchmarking against established planewave and real-space codes. Using this implementation, we study the band structure properties of X nanotubes and the bending properties of Xene sheets. Specifically, we first show that zigzag and armchair X nanotubes of radii from $1$ to $5$ nm are semiconducting, other than the armchair and zigzag type III carbon variants, for which we find a vanishingly small bandgap, indicative of metallic behavior. In particular, we find an inverse quadratic dependence of the bandgap with respect to radius for armchair and zigzag type III carbon nanotubes, other than which, there is an inverse linear dependence. Next, using the connection between cyclic symmetry and uniform bending deformations, we calculate the bending moduli of Xene sheets in both zigzag and armchair directions, while considering radii of curvature up to $5$ nm. For all sheets, we observe a Kirchhoff-Love type bending behavior, with graphene and stanene possessing the largest and smallest moduli, respectively. In addition, apart from graphene, there is significant anisotropy in the bending moduli, with larger values along the armchair direction. Finally, we demonstrate that the proposed method has very good parallel scaling and is highly efficient, enabling ab initio simulations of unprecedented size for systems with sufficiently high degree of cyclic symmetry, e.g., we simulate a silicon nanotube of radius $\sim 1\mu$m within $53$ minutes on $353$ processors.

The remainder of this manuscript is organized as follows. In Section II, we summarize the underlying real-space formulation of DFT that will be adopted here. In this framework, we develop a symmetry-adapted formulation of DFT for cylindrical geometries, as described in Section III. We then describe its parallel implementation in Section IV. Next, we demonstrate the accuracy and efficiency of the proposed formulation and implementation in Section V, a section in which we also study the electronic properties of X nanotubes and the bending properties of Xene sheets. Finally, we provide concluding remarks in Section VI.

II Real-space formulation of DFT

We first provide some mathematical background on Kohn-Sham DFT, while adopting a real-space formalism that has been shown to be both accurate and efficient Ghosh and Suryanarayana (2017a, b). Neglecting spin and utilizing the pseudopotential approximation, the electronic free energy of a system can be written as (Kohn and Sham, 1965; Mermin, 1965):

[TABLE]

where $\Psi$ is the collection of Kohn-Sham orbitals associated with the system, ${\mathbf{g}}$ is the corresponding collection of orbital occupation numbers, X is the collection of atomic positions, and $\rho$ is the electron density. The electron density can itself be expressed in terms of the orbitals and their occupations as:

[TABLE]

where $\bar{N}_{s}$ denotes the total number of orbitals. Above and in what follows, the generic index $n$ is used to label the orbitals and the corresponding orbital occupations (i.e., $\psi_{n}\in\Psi$ is a typical orbital and its occupation is $g_{n}\in{\mathbf{g}}$), while the generic index $J$ will provide labels for the nuclei (i.e., $\textbf{x}_{J}\in\textbf{X}$ is the position of the $J^{\text{th}}$ nucleus). In addition, the subscript $J$ will be used to explicitly indicate the dependence of a function (or a collection) on the species of the atom located at $\textbf{x}_{J}$.

The various terms arising in Eq. 1 can be interpreted as follows. The first term models the kinetic energy of the system of electrons, and it can be written as:

[TABLE]

The second term represents the exchange-correlation energy, for which many models exist, including the Local Density Approximation (LDA)(Kohn and Sham, 1965) and the Generalized Gradient Approximation (GGA) (Perdew et al., 1996). In this work, we employ the LDA:

[TABLE]

The third term accounts for the contribution from the nonlocal part of the pseudopotentials, which takes the following form within the Kleinman-Bylander representation Kleinman and Bylander (1982):

[TABLE]

Here, $\mathcal{A}_{J}$ denotes the collection of projectors associated with the atom at $\textbf{x}_{J}$, and $\chi_{J;p}$ are the nonlocal projection functions, with $\gamma_{J;p}$ representing the corresponding normalization constants. The fourth term represents the total electrostatic interaction energy, for which we use a local formulation of the electrostatics (Pask and Sterne, 2005b; Suryanarayana and Phanish, 2014; Ghosh and Suryanarayana, 2016) suitable for real-space DFT:

[TABLE]

Above, $\phi$ is the electrostatic potential, $b$ is the total pseudocharge density of the nuclei, and $\bar{E}_{sc}$ is the sum of the self energy and the repulsive energy corrections associated with the pseudocharges. Finally, the last term accounts for the electronic entropy arising from fractionally occupied electronic states at a given electronic temperature $T$:

[TABLE]

where $k_{B}$ is the Boltzmann constant.

The computation of the electronic ground state corresponding to the given set of (fixed) atomic positions X is given by the variational problem:

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta function and $\bar{N}_{e}$ is the total number of electrons. Note that the constraint of orthonormality on the orbitals is a consequence of the Pauli exclusion principle and the constraint on the occupations arises from the constancy of the number of electrons due to the Aufbau principle Ciarlet et al. (2003). In practice, the electronic ground state is often determined by seeking stationary states corresponding to solutions of the Euler-Lagrange equations:

[TABLE]

where $V_{xc}=\displaystyle\frac{\delta\bar{E}_{xc}}{\delta\rho}$ is the exchange-correlation potential, the electrostatic potential $\phi$ is the solution of the Poisson equation:

[TABLE]

the occupations $g_{n}$ are given by the Fermi-Dirac function:

[TABLE]

and $\bar{\mathcal{V}}_{nl}$—projection operator associated with the non-local part of the pseudopotential—acts on any given function $f({\mathbf{r}})$ as:

[TABLE]

Note that while solving Eqs. 9 and 10, the boundary conditions prescribed on the orbitals $\psi_{n}$ and the electrostatic potential $\phi$ are that they decay to zero at infinity (Hoffmann-Ostenhof et al., 1980; Ahlrichs et al., 1981; Suryanarayana et al., 2010). Here and henceforth, we assume that the systems are charge neutral.

Once the electronic ground-state has been determined, the free energy can be calculated using either Eq. 1 or the Harris-Foulkes Harris (1985); Foulkes and Haydock (1989) type functional :

[TABLE]

while the Hellmann-Feynman atomic forces required for geometry optimization and molecular dynamics can be written as:

[TABLE]

where $b_{J}$ represents the pseudocharge of the $J^{th}$ nucleus and $\text{Re}[.]$ denotes the real part of the bracketed expression. The first term is the local component of the force, the second term $\bar{\mathbf{f}}_{sc,J}(\textbf{X})=\frac{\partial\bar{E}_{sc}(\mathbf{X})}{\partial{\mathbf{x}}_{J}}$ corrects for overlapping pseudocharge densities, and the final term is the nonlocal component of the force.

III Symmetry-adapted real-space DFT on a cylinder

We now present a real-space formulation of DFT for cylindrical geometries that is able to explicitly incorporate the natural symmetries commonly arising in such systems. In order to achieve this, the Kohn-Sham equations posed on all of space in the previous section will be appropriately modified and augmented with suitable boundary conditions, as detailed below.

III.1 System specification: domain, atomic configuration and symmetries

Let $(\textbf{e}_{x},\textbf{e}_{y},\textbf{e}_{z})$ denote the canonical Cartesian coordinate axes and let $(r,\theta,z)$ represent the cylindrical coordinates associated with a point having Cartesian coordinates $(x,y,z)$. Consider an annular cylindrical region:

[TABLE]

with axis along $\textbf{e}_{z}$. The boundary of $\Omega$ can be decomposed as

[TABLE]

where $\partial R_{in}$ and $\partial R_{out}$ denote the surfaces $r=R_{in}$ and $r=R_{out}$, respectively; $\partial\vartheta_{0}$ and $\partial\vartheta_{\Theta}$ denote the surfaces $\theta=0$ and $\theta=\Theta$, respectively; and $\partial\mathcal{Z}_{0}$ and $\partial{\mathcal{Z}}_{H}$ denote the surfaces $z=0$ and $z=H$, respectively. Let this region $\Omega$ contain $N$ atoms positioned at ${\mathbf{R}}_{1},{\mathbf{R}}_{2},\ldots,{\mathbf{R}}_{N}$, the collection of which will henceforth be denoted by ${\mathbf{R}}$. A schematic illustrating this system is as shown in Fig. 1.

The system under consideration is allowed to have transational symmetry along the axis $\textbf{e}_{z}$ and/or rotational symmetry in the angular direction, i.e., cyclic symmetry about axis $\textbf{e}_{z}$. In the former case, the height of the domain $H$ must be commensurate with the translational periodicity of the system along $\textbf{e}_{z}$. In the latter, the maximum polar angle $\Theta$ associated with the domain must be of the form $\Theta=2\pi/\mathfrak{N}$, with $\mathfrak{N}$ being a natural number. In such cases, the region $\Omega$ is a fundamental domain or the unit cell for the symmetries involved, which plays a role similar to that of the unit cell employed in standard periodic calculations. In the current context, whenever periodic symmetries are present, the system under study is quasi-1D in nature. The symmetry group in this case is generated by translations along $\textbf{e}_{z}$, and is of the form:

[TABLE]

When cyclic symmetries are present, the symmetry group is generated by rotations with the common axis $\textbf{e}_{z}$, and it can be written as the set of matrices

[TABLE]

The symmetry groups associated with the (tubular) physical systems studied in this work are a combination of cyclic and translational symetries, i.e., the group in such cases is expressible as the direct product of the groups $\mathfrak{T}$ and $\mathfrak{C}$. Thus, the group is identifiable as the set of isometries (i.e., rigid body motions):

[TABLE]

and it can be indexed by pairs of numbers $(\zeta,\mu)$ with $\zeta=0,1,2,\ldots,\mathfrak{N}-1$ and $\mu\in\mathbb{Z}$. Specifically, the group element associated with the pair $(\zeta,\mu)$ is the isometry ${\Upsilon}_{\zeta,\mu}=(\mathfrak{R}^{\zeta}|\textbf{t}_{\mu})$, whose action on a point in space (denoted henceforth as $\circ$) rotates it by $\mathfrak{R}^{\zeta}$ about $\textbf{e}_{z}$, while also simultaneously translating it by $\mu H$ along $\textbf{e}_{z}$. The global tubular structure X that is effectively being simulated can be generated as the image of the points ${\mathbf{R}}_{1},{\mathbf{R}}_{2},\ldots,{\mathbf{R}}_{N}$ under the action of the isometries in the group $\mathcal{G}$, i.e.,

[TABLE]

Correspondingly, the global simulation domain $\mathcal{C}$, that encases all of the points in X, can be generated as the image of the fundamental domain $\Omega$ under the group $\mathcal{G}$. We will use the notation $\mathcal{G}\circ{\mathbf{R}}_{J}$ to denote the orbit of the point ${\mathbf{R}}_{j}$ under the group, i.e.,

[TABLE]

With this notation, the collection of points X is expressible as:

[TABLE]

Indeed, in the context of DFT, X is nothing but the collection of atomic positions in ${\mathbb{R}}^{3}$.

III.2 Formulation: Symmetry-adapted real-space DFT

A basic consequence of the presence of physical symmetries in a system—specifically, the atomic positions of the structure being describable as the orbit of a discrete group of isometries—is that, under some generally applicable hypotheses (Banerjee, 2013; Banerjee and Suryanarayana, 2016), the electron density for such a system is invariant under the symmetry group and further, the Kohn-Sham Hamiltonian for the system commutes with the symmetry operations of the group. Results from group representation theory dictate therefore, that the eigenstates of this operator can be characterized through the irreducible representations of the symmetry group, and that the eigenstates transform as the irreducible representations under the action of the group(Banerjee, 2013; McWeeny, 2002; Hamermesh, 2012). The relevant symmetry group $\mathcal{G}$ in the present context is Abelian, therefore, its complex irreducible representations are one dimensional (Folland, 1994; Barut and Raczka, 1986). These one-dimensional irreducible representations, or the so called complex characters of $\mathcal{G}$ are complex valued functions of the group: identifying the group element $\Upsilon_{\zeta,\mu}\in\mathcal{G}$ in terms of the pair $(\zeta,\mu)\in\{0,1,2,\ldots,\mathfrak{N}-1\}\times\mathbb{Z}$, the set of characters of $\mathcal{G}$ can be written as:

[TABLE]

The variables ${\nu}$ and ${\eta}$ serve to label the complex characters of $\mathcal{G}$, and consequently, they also label the associated eigenstates of the Kohn-Sham Hamiltonian. In what follows, we explicitly indicate this labeling for eigenvalues, eigenvectors, and occupations as $\lambda_{n}(\nu,\eta),\psi_{n}({\mathbf{r}},\nu,\eta)$, and $g_{n}(\nu,\eta)$ respectively.

A number of consequences are associated with, or result from the above observations, as we now discuss. First, due to the orthogonality relations obeyed by the characters(Folland, 1994; Barut and Raczka, 1986), the collections of eigenstates associated with distinct characters are mutually orthogonal. Hence, by the use of a symmetry adapted basis (McWeeny, 2002), the Hamiltonian $\bar{\mathcal{H}}$ can be block-diagonalized(Banerjee, 2013) and the eigenvalue problems associated with distinct characters (i.e., distinct values of $(\eta,\nu)$) can be solved independently of one another. Second, the boundary conditions on the orbitals, and the electrostatic potential, that have to be applied on certain surfaces of the computational domain $\Omega$ can be readily deduced based on transformation properties of the characters. Third, any quantity that involves contributions from all eigenstates that appear in the problem, has to include contributions from each of the elements of $\widehat{\mathcal{G}}$ – this can be achieved by integrating relevant eigenstate-dependent quantities against a suitable integration measure over $\widehat{\mathcal{G}}$. As an example, consider the electron density: if each of the diagonal blocks of the symmetry-adapted Hamiltonian contributes $N_{s}$ electronic states, Eq. 2 can be rewritten as:

[TABLE]

Here, the sum $\displaystyle\frac{1}{\mathfrak{N}}\sum_{\nu=0}^{\mathfrak{N}-1}$ is associated with integrating against $\nu$, and $\displaystyle\fint$, which signifies the average over the interval $\displaystyle\left[-\frac{\pi}{H},\frac{\pi}{H}\right]$, accumulates contributions in $\eta$. We now exploit these consequences to reduce the Kohn-Sham problem described in Section II to the fundamental domain $\Omega$.

III.2.1 Boundary conditions

Orbitals

The eigenfunctions of the Hamiltonian transform in accordance with the irreducible representations, whereby the action of an arbitrary group element $\Upsilon_{\zeta,\mu}\in\mathcal{G}$ on an orbital associated with the character $(\nu,\eta)$ can be written as:

[TABLE]

or equivalently

[TABLE]

These can be identified as versions of the Bloch-theorem(Bloch, 1929; Odeh and Keller, 1964) associated with the symmetry group $\mathcal{G}$. Therefore, we arrive at the following boundary conditions for the orbitals on the surfaces $\partial\vartheta_{0}{\,\textstyle\bigcup\,}\partial\vartheta_{\Theta}$ and $\partial\mathcal{Z}_{0}{\,\textstyle\bigcup\,}\partial{\mathcal{Z}}_{H}$, respectively:

[TABLE]

For the surfaces $\partial R_{in}$ and $\partial R_{out}$, we assume that the atoms within $\Omega$ are sufficiently far from these surfaces, allowing the decay of the electron density along the radial direction to come into effect, i.e.,

[TABLE]

Electrostatic potential

The electron density $\rho$ is group invariant, i.e., it is transforms under the group as functions associated with the characters $(\nu,\eta)=(0,0)$. It can be easily shown that the total pseudocharge density $b$ inherits this symmetry as well (Banerjee and Suryanarayana, 2016), since the atomic positions are expressible as the orbit of the group and the individual pseudocharges are spherically symmetric. Therefore, it follows from Eq. 10 that $\phi$ is group invariant, which implies that on $\partial\vartheta_{0}{\,\textstyle\bigcup\,}\partial\vartheta_{\Theta}$ and $\partial\mathcal{Z}_{0}{\,\textstyle\bigcup\,}\partial{\mathcal{Z}}_{H}$, respectively, we have:

[TABLE]

For the surfaces $\partial R_{in}$ and $\partial R_{out}$, the boundary conditions can be determined by using the integral form (Evans, 1998) of the solution to Eq. 10:

[TABLE]

which can then be evaluated using Ewald summation(Miller and Rouet, 2010; Langridge et al., 2001) or multipole expansion(Han et al., 2008; Ghosh and Suryanarayana, 2017a) techniques.

III.2.2 Energy, Kohn-Sham equations, and atomic forces

In the discussion that follows, we denote the collection of character dependent electronic states by $\Psi(\widehat{\mathcal{G}})$ and the corresponding collection of electronic occupations by ${\mathbf{g}}(\widehat{\mathcal{G}})$. If, as before, each diagonal block of the Kohn-Sham Hamiltonian contributes $N_{s}$ states, we have:

[TABLE]

Energy functional

The presence of symmetries make the global system extended in nature, i.e., the atoms within the simulation domain $\Omega$ represent only a portion of the global infinite structure. Therefore, the energies have to be interpreted in a per fundamental domain sense. Moreover, although the computation is confined to the fundamental domain, the various terms in Eq. 1 have to be suitably modified to account for (i) the effect of atoms which belong to the global structure but lie outside the fundamental domain, and (ii) the (possibly infinite) multiplicities of electronic states arising due to the effects of symmetry and the extended nature of the system. Keeping these in mind, the symmetry-adapted electronic free energy per fundamental domain can be written as

[TABLE]

whose terms are now described in detail.

The first term in the energy functional (Eq. 33) is the electronic kinetic energy per fundamental domain, and like the electron density, it includes contribution from $N_{s}$ electronic states from each diagonal block of the Hamiltonian. As a result, Eq. 3 is modified to:

[TABLE]

The second term in the energy functional (Eq. 33) is the exchange-correlation energy per fundamental domain, and since the electron density obeys the symmetry of the structure, Eq. 4 reduces to:

[TABLE]

Note that even though we are focusing on LDA in this work, an analogous expression involving the gradient of the electron density is applicable for semilocal exchange-correlation functionals such as the GGA.

The third term in the energy functional (Eq. 33) is the nonlocal pseudopotential energy per fundamental domain. To obtain this term from Eq. 5, we include contribution from $N_{s}$ electronic states from each diagonal block of the Hamiltonian to arrive at:

[TABLE]

Here, we have accumulated the contribution of the projectors centered on the atoms within the fundamental domain, as applied to all the electronic states in the system. Note that since the atom centered projectors can have support extending beyond the fundamental domain, their overlaps with the orbitals need to be carried out over the global simulation domain $\mathcal{C}$. Eq. 36 can now be rewritten as follows

[TABLE]

where the function:

[TABLE]

The first equality in Eq. 37 has been obtained by expressing the the global simulation domain $\mathcal{C}$ as the image of the fundamental domain $\Omega$ under the action of the group $\mathcal{G}$. The second equality is obtained by employing the change of variables ${\mathbf{r}}=\Upsilon_{\zeta,\mu}\circ\textbf{y}$ and then using Eq. 26. To obtain the third equality, it should be noted that the atomic position ${\mathbf{R}}_{J}$ can be expressed as $\Upsilon_{\zeta,\mu}\circ{\mathbf{R}}_{J^{\prime}}$, with ${\mathbf{R}}_{J^{\prime}}$ denoting an image of ${\mathbf{R}}_{J}$ that (in general) lies away from the fundamental domain. Correspondingly, the quantities $\displaystyle\frac{-2\pi\zeta}{\mathfrak{N}}$ and $-\mu H$ that appear in the exponential can be expressed in terms of the differences in polar coordinates $(r_{J},\theta_{J},z_{J})$ and $(r_{J^{\prime}},\theta_{J^{\prime}},z_{J^{\prime}})$ of the atomic positions ${\mathbf{R}}_{J}$ and ${\mathbf{R}}_{J^{\prime}}$. In particular, we have used (Banerjee and Suryanaryana., 2019) the fact that the projectors $\chi_{J,p}$ have a mathematical form that is similar to atomic orbitals, i.e, a spherical harmonic multiplied by a compactly supported radial function.

The fourth term in the energy functional (Eq. 33) is the electrostatic interaction energy per fundamental domain. Since all quantities involved obey the symmetry of the structure, it is sufficient to work with their restriction to the fundamental domain. Therefore, Eq. 6 reduces to:

[TABLE]

where the pseudocharge density $b({\mathbf{r}},\textbf{R},\mathcal{G})$ admits the decomposition:

[TABLE]

Note that the term $\bar{E}_{sc}$ can similarly be reduced to $E_{sc}$ by restricting the integrals over ${\mathbb{R}}^{3}$ to $\Omega$, not described here for brevity.

The fifth and final term in the energy functional (Eq. 33) represents the electronic entropy contribution to the free energy per fundamental domain. To obtain this term from Eq. 7, we include contribution from the occupations of $N_{s}$ electronic states from each diagonal block of the Hamiltonian to arrive at:

[TABLE]

Variational problem

The symmetry-adapted variational problem for determining the electronic ground state takes the form:

[TABLE]

where the symmetry-adapted orthogonality constraints follow from the orthogonality of the orbitals associated with the distinct characters of $\mathcal{G}$. The constraint on the number of electrons $N_{e}$ in the fundamental domain is obtained by including the contribution from the occupations of $N_{s}$ electronic states from each diagonal block of the Hamiltonian.

Kohn-Sham equations

On taking variations of the constrained minimization problem presented above, we arrive at the following symmetry-adapted Kohn-Sham equations on the fundamental domain $\Omega$:

[TABLE]

subject to the boundary conditions given by Eqs. 27-29. In the above equation, $V_{xc}=\displaystyle\frac{\delta E_{xc}}{\delta\rho}$ is the exchange-correlation potential, $\phi$ is the solution of the following Poisson problem on $\Omega$:

[TABLE]

subject to the boundary conditions given by Eqs. 30-32. The symmetry-adapted nonlocal operator $\mathcal{V}_{nl}$ acts on a function $f$ defined over the fundamental domain as:

[TABLE]

In the above equations and those that follow, the electron density is calculated using Eq. 24, and the occupations are determined through the Fermi-Dirac function, i.e.,

[TABLE]

where the Fermi level $\lambda_{F}$ is determined by satisfying the constraint on the number of electrons in the fundamental domain (Eq. III.2.2).

Harris-Foulkes functional

Once the electronic ground state is determined by solving the above equations self-consistently, the free energy per fundamental domain can be computed through Eq. 33, or alternately through the symmetry-adapted Harris-Foulkes functional:

[TABLE]

Atomic forces

The symmetry-adapted Hellmann-Feynman atomic forces in the Cartesian coordinate system take the form:

[TABLE]

where $\nabla$ denotes the Cartesian gradient operator and $\mathbf{f}_{sc,J}({\mathbf{R}},\mathcal{G})=-\frac{\partial E_{sc}({\mathbf{R}},\mathcal{G})}{\partial{\mathbf{R}}_{J}}$. Note that in deriving the nonlocal component of the force, we have transferred the derivative on the projectors (with respect to atomic position) onto the orbitals (with respect to atomic position). This is because the orbitals are typically smoother than the projector functions, therefore resulting in substantially more accurate atomic forces (Hirose et al., 2005; Ghosh and Suryanarayana, 2017a, b).

III.2.3 Time reversal symmetry

In the absence of magnetic fields, the symmetry-adapted formulation presented above can be further reduced by employing time-reversal symmetry(Martin, 2004; Banerjee and Elliott., 2019). Specifically, for $\nu\in\{1,2,\ldots,\mathfrak{N}-1\}$, we have the relations:

[TABLE]

and for $\nu=0$:

[TABLE]

As a result, the character space $(\nu,\eta)$ that needs to be considered for the Kohn-Sham equations as well as the calculation of the free energy and atomic forces is essentially halved.

III.2.4 Computational cost reduction due to cyclic symmetry adaptation

The proposed formulation is able to achieve significant reduction in the computational cost of DFT simulations for systems possessing cyclic symmetry. Consider the nonlinear eigenvalue problem, which forms the dominant part of Kohn-Sham calculations and scales cubically with the problem size asymptotically. Since the problems associated with different $\nu$ are independent, the cyclic symmetry adaptation translates to a factor of $\mathcal{O}(\mathfrak{N}^{2})$ reduction in computational cost. This can result in tremendous savings, particularly for systems where $\mathfrak{N}$ is large, such as those studied in this work. Reductions in cost also extend to the solution of the Poisson equation as well as the calculation of the energy and forces, where the factor is the more modest $\mathcal{O}(\mathfrak{N})$. Note that even further reductions are achieved in practice while solving for the electronic ground state using the self-consistent field (SCF) method, due to the notable reduction in symmetry breaking and charge sloshing type instabilities Banerjee and Suryanarayana (2016). In addition, the implementations are more amenable to scalable parallel computations, due to the significantly fewer global communications. Indeed, such communications are significantly reduced because the orbitals associated with distinct characters are automatically orthogonal to each other.

IV Numerical Implementation

We implement the proposed formulation within the real-space finite-difference DFT code SPARC Ghosh and Suryanarayana (2017a, b). Due to the nature of SPARC’s implementation, we employ the Bloch-type ansatz for convenience:

[TABLE]

where $u_{n}({\mathbf{r}},\nu,\eta)$ is group invariant, i.e., for any $\Upsilon_{\zeta,\mu}\in\mathcal{G}$,

[TABLE]

As a result, $u_{n}$ replaces $\psi_{n}$ as the primary unknown functions that are being solved for in the symmetry-adapted Kohn-Sham equations, the modified form of which can easily be derived using the above ansatz, and are therefore not reproduced here for the sake of brevity.

The geometry of $\Omega$ motivates the use of cylindrical polar coordinates, wherein the Laplacian and the Cartesian gradient operator take the form:

[TABLE]

The fundamental domain $\Omega$ is discretized using a finite-difference grid with spacing $h_{r}$, $h_{\theta}$, and $h_{z}$ along the $r$, $\theta$, and $z$ directions, respectively. This implies that $R_{2}-R_{1}=n_{r}h_{r}$, $\Theta=n_{\theta}h_{\theta}$ and $H=n_{z}h_{z}$, for natural numbers $n_{r},n_{\theta}$ and $n_{z}$. Each finite-difference node is indexed using a triplet of the form $(i,j,k)$, with $i=1,2,\ldots,n_{r}$, $j=1,2,\ldots,n_{\theta}$, and $k=1,2,\ldots,n_{z}$. Using the central finite-difference approximation, we approximate the partial first derivatives as:

[TABLE]

and similarly the partial second derivatives as:

[TABLE]

with $f^{(i,j,k)}$ representing the value of the function $f$ at the node $(i,j,k)$. Denoting $s\in\{r,\theta,z\}$, the weights that appear in the above expressions can be written as (Mazziotti, 1999; Suryanarayana and Phanish, 2014):

[TABLE]

Due to the curvilinear nature of the underlying coordinate system, the Laplacian and Hamiltonian matrices resulting from the above discretization scheme are non-Hermitian, even though the infinite-dimensional operators from which they arise are Hermitian(Banerjee and Suryanarayana, 2016; Gygi and Galli, 1995). For a fixed finite-difference order however, as the discretization is refined, the discrete Laplacian and Hamiltonian matrices approach Hermitian matrices and the eigenvalues of these matrices turn out to be either real, or they have vanishingly small imaginary partsGygi and Galli (1995); Banerjee and Suryanarayana (2016). Hence this issue does not negatively impact the physical results obtained or their implications.

We approximate the integrals over the fundamental domain by employing the following quadrature rule:

[TABLE]

with $r_{i}$ denoting the radial coordinate of the finite-difference node indexed by $(i,j,k)$. We enforce periodic boundary conditions by mapping any index that does not correspond to a node in the finite-difference grid to its periodic image within $\Omega$. We enforce zero Dirichlet boundary conditions by setting $f^{(i,j,k)}=0$ for any index that does not correspond to a node in the finite-difference grid. Following the strategy proposed previouslySuryanarayana et al. (2013); Ghosh and Suryanarayana (2017a), we use the discrete Laplacian to directly compute the pseudocharges from the local parts of the pseudopotentials, while assigning them to the grid. Due to the presence of translational symmetry along $\textbf{e}_{z}$, evaluating integrals such as Eq. 24 requires us to discretize the domain of the variable $\eta$, allowing for a discrete representation of the complex characters associated with the periodic symmetry. Accordingly, we utilize the Monkhorst-Pack Monkhorst and Pack (1976) grid for sampling the interval $[-\frac{\pi}{H},\frac{\pi}{H}]$ and approximate the averaged integral of any function over the interval as:

[TABLE]

Here $\eta_{b}$ and $w_{b}$ denote the integration nodes and weights, respectively. The total number of discretized characters in the computation (i.e., “k-points” in the language of periodic DFT calculations) is denoted as $N_{K}\approx\frac{1}{2}(\mathfrak{N}\times N_{\eta})$, where time-reversal symmetry has been used to reduce the total number by a factor of two approximately.

We use the Chebyshev polynomial filtered subspace iteration (CheFSI) technique Zhou et al. (2006a, b) in conjunction with potential mixing for computing the electronic ground state corresponding to Eq. 43. Within the CheFSI method, we employ Arnoldi iterations(Saad, 2003) for calculating the extremal eigenvalues of the Hamiltonian, and LAPACKAnderson et al. (1999) for solving the projected subspace eigenproblem. We solve the linear system corresponding to the Poisson problem (Eq. 44) using the block-Jacobi preconditioned Golub and Van Loan (2012) Generalized minimal residual method (GMRES) Saad and Schultz (1986). We calculate the Fermi energy using Brent’s method Press (2007), and use Periodic Pulay extrapolation Banerjee et al. (2016) for accelerating convergence of the SCF iterations. We calculate the ground state free energy using the symmetry-adapted Harris-Foulkes type functional (Eq. 47). If and when required, we employ the FIRE algorithm Bitzek et al. (2006) for performing structural relaxations.

The real-space discretization naturally lends itself to parallelization via domain decomposition. In addition, the eigenvalue problems corresponding to distinct characters can be solved independently. Therefore, we employ two levels of parallelization: the first being over the different values of the characters and the second being over the spatial domain. The former is achieved by uniformly distributing the list of $N_{K}$ eigenvalue problems among $N^{proc}_{K}$ processors ($N^{proc}_{K}\leq N_{K}$). The latter is achieved by partitioning the fundamental domain $\Omega$ among $N^{proc}_{D}$ processors as:

[TABLE]

and assigning the portion of the calculation associated with the partition $\Omega_{p}$ to the $p^{th}$ processor. The total number of processors employed in this two level parallelization scheme is therefore $N^{proc}=N^{proc}_{D}\times N^{proc}_{K}$. We use the Portable, Extensible Toolkit for Scientific computations (PETSc)(Balay et al., 2018) suite of data structures and routines, in conjunction with the Message Passsing Interface (MPI)(Gropp et al., 1999) for implementation and parallelization of our computational routines.

V Results and discussion

The main objects of study in this work are single walled nanotubes of carbon, silicon, germanium, and tin, collectively referred to here as X (X=C,Si,Ge,Sn) nanotubes. These 1D nanostructures are formed by rolling their 2D sheet counterparts: graphene, silicene, germanene and stanene, collectively referred to here as Xenes. Depending on whether the direction of rolling is armchair or zigzag, the nanotubes can be classified as armchair or zigzag, respectively. Since zigzag carbon nanotubes have distinct electronic properties based on their radius (Saito et al., 1998), we further classify the zigzag X nanotubes as type I, II, or III, depending on whether $\text{mod}(\mathfrak{N},3)=1,2$ or [math]. Both X nanotubes and Xene sheets are known to demonstrate unusual and fascinating material properties Martel et al. (1998); Javey et al. (2003); Popov (2004); Gong et al. (2009); Park et al. (2009); Wu et al. (2012); Park et al. (2011); Li et al. (2011); Zhao et al. (2006); Xu et al. (2013); Bhimanapati et al. (2015); Butler et al. (2013); Naguib et al. (2014); Fiori et al. (2014); Koppens et al. (2014), motivating their choice here as well as in a number of previous electronic structure studies Blase et al. (1994); Spataru et al. (2004); Yang and Han (2000); Fagan et al. (2000); Benedict et al. (1995); Zhang et al. (2003); Yang and Ni (2005); Giovannetti et al. (2008); Vogt et al. (2012); Dávila et al. (2014); Zhu et al. (2015).

The calculations here utilize the LDA (Kohn and Sham, 1965) to model the exchange-correlation functional, with the Perdew-Wang parametrization Perdew and Wang (1992) of the correlation energy as calculated by Ceperley and Alder Ceperley and Alder (1980); smearing of $k_{B}T=0.001$ Ha, treated here as a numerical parameter to aid SCF convergence rather than the actual temperature; and Troullier-Martins norm conserving pseudopotentials(Troullier and Martins, 1991). With these choices, the ground state interatomic distance ($a$), and out of plane buckling distance ($\delta$) for the planar Xene sheets are as reported in Table 1. The agreement of these values with the literature is generally quite good, thus giving us confidence in the quality of the simulations.

The combined use of cyclic and periodic symmetries allows X nanotubes to be represented by just 4 atoms within the fundamental domain, i.e., the fundamental domain corresponds to the rolling of the 4-atom orthogonal unit cell in the Xene sheet, as shown in Fig. 2. The angle formed by the fundamental domain $\Theta=2\pi/\mathfrak{N}$, where $\mathfrak{N}$ depends on the the radius of the nanotube and the interatomic distance in the flat sheets. Specifically, in the absence of relaxation effects 111X (X=C, Si, Ge, Sn) nanotubes formed by rolling the corresponding Xene sheets are not necessarily at the structural ground state Sánchez-Portal et al. (1999). Indeed, for small radii nanotubes, the atoms can experience atomic forces as large as $0.01$ Ha/Bohr. However, for the large radii nanotubes studied in this work, the maximum component of the atomic force in the initial configuration is less than $3\times 10^{-3}$ Ha/Bohr, with the value becoming smaller as the radius gets larger. Therefore, the adopted procedure provides a very good guess for the atomic positions. Moreover, since there is no noticeable change in the results between this configuration and the structural ground state, we do not perform any structural optimization steps in this work. , $\mathfrak{N}=\pi/\sin^{-1}\left(\frac{L}{2R}\right)$, where $L=3a$ and $\sqrt[]{3}a$ for armchair and zigzag nanotubes, respectively. The corresponding heights of the fundamental domain are $H=\sqrt[]{3}a$ and $3a$, respectively. The radii $R_{in}$ and $R_{out}$ are chosen such that all atoms are at least $11$ Bohr away from the boundaries in the radial direction, so as to allow sufficient decay of the electron density and orbitals.

We discretize the governing equations using a twelfth-order accurate finite-difference discretization. Subsequent to a convergence analysis with respect to mesh size h=\displaystyle\max\left\{h_{r},\bigg{(}\frac{R_{in}+R_{out}}{2}\bigg{)}h_{\theta},h_{z}\right\} (similar to Section V.1) as well as an increasingly finer sampling of the values of $\eta$, the $h$ and $N_{\eta}$ listed in Table 2 are chosen for the X nanotube simulations in Sections V.2 and V.3, from which the physical properties of interest, namely bandgap variation with radius of X nanotubes and bending moduli of Xene sheets, are calculated. The chosen parameters ensure that the energy and atomic forces are converged to within $10^{-5}$ Ha/atom and $10^{-5}$ Ha/Bohr, respectively. Note that such high precision—significantly more stringent than that typically employed in DFT calculations—is essential to capture the small bandgap and energy variations that occur with respect to the radius.

V.1 Convergence and accuracy

To assess the accuracy of the proposed formulation and implementation, we first verify the convergence of the energy as well as the atomic forces with respect to spatial discretization, i.e., mesh size $h$. As representative systems, we choose zigzag X (X=C,Si,Ge,Sn) nanotubes with radii $0.90$, $0.98$, $0.99$, and $1.13$ nm, corresponding to $\mathfrak{N}=23$, $16$, $16$, and $16$, respectively. Note that the radii of these tubes is sufficiently small for the atoms to experience significant atomic forces. Without loss of generality, only the $\eta=0$ point is included for this numerical test (equivalent to a $\Gamma$-point calculation in traditional DFT). It is clear from the results in Fig. 3 that there is systematic convergence of both the energy and atomic forces to reference values obtained for $h=0.1$ Bohr. On fitting the data, we find average convergence rates in the energy and atomic forces of $5.5$ and $7.8$, respectively, comparable to those obtained by the analogous real-space formalism for affine coordinate systems Ghosh and Suryanarayana (2017a, b); Sharma and Suryanarayana (2018b).

In order to further verify the accuracy of the proposed method, we compare the results at $h=0.1$ Bohr with highly converged values obtained by the established planewave code ABINIT (Gonze et al., 2009, 2002). We find that there is agreement to within $6\times 10^{-5}$ Ha/atom and $1\times 10^{-4}$ Ha/Bohr in the energy and atomic forces, respectively. Indeed, even better agreement would have been possible, but for the significant challenge in converging ABINIT to finer levels of accuracy. This is mainly due to the stagnation in the results with respect to vacuum Ghosh and Suryanarayana (2017b), likely due to the inaccurate electrostatics resulting from the requirement of periodic boundary conditions. To confirm this, we have also compared with the real-space DFT code SPARC Ghosh and Suryanarayana (2017a, b), which is not only more efficient, but also does not suffer from the aforementioned stagnation. We have found that there is indeed better agreement with SPARC, with energy and atomic forces differing by not more than $2\times 10^{-6}$ Ha/atom and $1\times 10^{-5}$ Ha/Bohr, respectively.

Even though not demonstrated here, we have verified that the proposed method inherits a number of attractive features of the underlying SPARC framework. In particular, there is exponential convergence in the properties of interest with respect to the amount of vacuum in the radial direction. In addition, the computed energy and atomic forces are consistent, which allows for accurate geometry optimization and molecular dynamics simulations. Finally, the eggbox effect resulting from breaking of the translational and cyclic symmetry of the system is negligible, particularly at the mesh sizes considered in this work. Since these features have been demonstrated and discussed in detail previously for SPARC Ghosh and Suryanarayana (2017a, b), we do not repeat them here for the sake of brevity.

V.2 Band structure of X (X=C, Si, Ge, Sn) nanotubes

We now use the proposed method to study the band structure of X (X=C, Si, Ge, Sn) nanotubes, with radii ranging from $R=1.3$ nm ($\mathfrak{N}=23$) to $R=4.8$ nm ($\mathfrak{N}=108$). Studies in literature have shown that larger radii carbon nanotubes are more prone to instability (e.g., flattening or collapse) when subject to hydrostatic stresses. Chopra et al. (1995); Gao et al. (1998); Elliott et al. (2004); Tangney et al. (2005); Tang et al. (2005); Gadagkar et al. (2006). There is however substantial disagreement in the theoretically/computationally predicted values, with critical radii ranging from $1$ to $3.5$ nm at atmospheric conditions. Concerns regarding the validity of these predictions remain, particularly with the synthesis of carbon nanotubes having radii as large as $6$ nm Cheung et al. (2002). The systems chosen here are motivated by the fact that large radii nanotubes are far less studied, particularly in the context of ab-initio calculations, where the associated computational cost is large. Moreover, such radii are required for accurately calculating the bending moduli of the Xene sheets, as done in Section V.3.

We start by computing the symmetry-adapted band structure data for the aforementioned nanotubes in the discrete $(\nu,\eta)$ space. This ability to calculate and plot the variation of the eigenvalues with respect to both the character labels $(\nu,\eta)$, rather than with respect to $\eta$ alone (as is typically done), is a distinctive feature of the method developed here. Such symmetry-adapted band structure diagrams have the advantage that they allow for significantly easier interpretation of the results, particularly for systems with complex band structure. See Fig. 4 for representative band structure diagrams of carbon armchair ($R=1.6$ nm, $\mathfrak{N}=23$) and tin zigzag type I ($R=1.5$ nm, $\mathfrak{N}=22$) nanotubes along specific line segments in $(\eta,\nu)$ space. Note that for the band structure diagrams at fixed $\eta$ (i.e., Figs. 4b and 4d), $\nu$ can only take integer values.

The calculated band structure data can be used to deduce whether the systems are metallic, semi-metallic or insulating. For insulating systems, the size of the bandgap can be determined by calculating the difference between the smallest eigenvalue above the Fermi level and largest eigenvalue below the Fermi level in all of $(\eta,\nu)$ space. We have found that all the X nanotubes are semiconducting. Specifically, there is a direct bandgap for the armchair nanotubes at $\left(\frac{\eta H}{2\pi},\nu\right)=\left(\frac{1}{3},0\right)$ (or equivalently $\left(-\frac{1}{3},0\right)$). In addition, apart from the zigzag type I carbon nanotube which has a direct bandgap at $\left(\frac{\eta H}{2\pi},\nu\right)=\left(0,\frac{\mathfrak{N}-1}{3}\right)$ (or equivalently $\left(0,\frac{2\mathfrak{N}+1}{3}\right)$), the other zigzag nanotubes have a direct bandgap at $\left(\frac{\eta H}{2\pi},\nu\right)=\left(0,\frac{\mathfrak{N}+2}{3}\right)$ (or equivalently $\left(0,\frac{2\mathfrak{N}-2}{3}\right)$), $\left(\frac{\eta H}{2\pi},\nu\right)=\left(0,\frac{\mathfrak{N}+1}{3}\right)$ (or equivalently $\left(0,\frac{2\mathfrak{N}-1}{3}\right)$), and $\left(\frac{\eta H}{2\pi},\nu\right)=\left(0,\frac{\mathfrak{N}}{3}\right)$ (or equivalently $\left(0,\frac{2\mathfrak{N}}{3}\right)$) for the type I, II, and III variants, respectively. 222The implication (if any) of the bandgap being at different values of $\nu$ is not evident at this point, and requires a further in-depth study.

The above results are generally in good agreement with those found in literature Hamada et al. (1992); Saito et al. (1992); Mintmire et al. (1993); Wang et al. (2017); Yang and Ni (2005), apart from a few discrepancies. First, the bandgap in armchair germanium and tin nanotubes has previously been found to be indirect Wang et al. (2017), whereas we predict a direct bandgap. This is likely due to the much larger nanotubes studied here. Second, while some electronic structure studies Wang et al. (2017), including the one here, have shown armchair silicon nanotubes to be semiconducting, others have found them to be metallic Fagan et al. (2000); Yang and Ni (2005). A possible source of the disagreement is that rather small planewave cutoffs have been used in cases where metallic behavior has been predicted. Finally, armchair carbon nanotubes are generally found to be metallic Hamada et al. (1992); Saito et al. (1992); Mintmire et al. (1993), however here we have obtained a nonzero (but vanishingly small) bandgap. We have verified that the above disagreements are not an artifact of the proposed formulation or the chosen pseudopotential/exchange-correlation functional. For example, consider the armchair carbon nanotube of radius $R=1.55$ nm. Highly accurate planewave calculations using ABINIT predicts a bandgap identical to the value computed here (i.e., 0.0136 eV). On changing the pseudopotential from Troullier-Martins to ONCV Hamann (2013), the bandgap remains with a nearly identical value of 0.0133 eV, and on further changing the exchange-correlation functional from LDA to GGA Perdew et al. (1996), the bandgap persists with a nearly identical value of 0.0130 eV. Since it is known from symmetry arguments that these nanotubes are metallic (Saito et al., 1998), the vanishingly small bandgaps observed in the simulations are possibly a consequence of symmetry breaking arising due to numerical artifacts. In any case, given the negligible bandgaps, the results here indicate that armchair nanotubes are metallic at ambient conditions White and Mintmire (2005), in agreement with previous work as well as experimental measurements Tans et al. (1997).

Next, we determine the variation of bandgap with nanotube radius $R$, the results of which are presented in Fig. 5. Anticipating an inverse power-law dependence, we compute the decay exponents through straight line fits and present the results so obtained in Table 3. It is clear that the nearly all nanotubes possess a close to inverse linear dependence with radius, the exceptions being armchair and zigzag type III nanotubes of carbon, which possess a close to inverse quadratic dependence. This atypical dependence in carbon nanotubes is consistent with results obtained from elaborately constructed tight binding models for graphene that are able to explicitly account for curvature effects (Ding et al., 2002). Note that these effects are automatically incorporated into our ab-initio simulations and therefore provide an elegant route to the fitting of the material parameters that appear in such tight binding models. The minor deviation of the computed decay exponents from $1$ or $2$ is perhaps indicative that the bandgap of such materials is better expressed by a relationship of the form $\frac{c_{1}}{R}+\frac{c_{2}}{R^{2}}$. However, these variations are quite small for the Xene tubes considered here, and therefore have been neglected. Overall, apart from armchair carbon nanotubes that have been previously considered as metallic, the scaling laws obtained here for the bandgap as a function of the radius are in good agreement with literature (Odom et al., 2000; Ding et al., 2002; Ouyang et al., 2002; Wang et al., 2017)

It is worth noting that the LDA exchange-correlation functional has limitations in the prediction of quantitatively accurate bandgaps (Sham and Schlüter, 1985; Hybertsen and Louie, 1985; Perdew and Levy, 1983; van Schilfgaarde et al., 2006). However, qualitative trends such as those discussed above are expected to be representative of the physical behavior. Indeed, the proposed formulation does not have any fundamental difficulty in dealing with semilocal and hybrid exchange-correlation functionals, which can be employed for making more quantitatively accurate predictions, making it a worthy subject for future work.

V.3 Bending moduli of Xene (X=C, Si, Ge, Sn) sheets

We now use the proposed method to calculate the bending moduli of Xene (X=C, Si, Ge, Sn) sheets along the armchair and zigzag directions. Specifically, we consider uniformly bent sheets with radii of curvature ranging from $R=1.3$ nm to $R=4.8$ nm. In order to significantly increase the efficiency of the simulations, we approximate a Xene sheet bent along the armchair or zigzag directions as armchair or zigzag nanotubes, respectively, with the radius of the nanotube chosen such that it matches the desired radius of curvature.Banerjee and Suryanarayana (2016) This strategy indeed neglects edge related effects, which are expected to play a relatively minor role in the current context. Such an approximation can be justified by appealing to Saint-Venant’s principle (Iesan, 2006) and the nearsightedness of matter (Prodan and Kohn, 2005) at the continuum and electronic structure scales, respectively.

Using the strategy described above, we calculate the bending energy333In the context of nanotubes, this is referred to as formation energy or strain energy. $\mathcal{E}_{bend}$ as a function of the radius of curvature $R$ along the armchair and zigzag directions, and plot the results so obtained in Fig. 6. The bending energy $\mathcal{E}_{bend}$ at a given radius of curvature $R$ is defined to be the difference in the free energy per fundamental domain between the nanotube (with radius $R$) and sheet configurations, normalized by the area of the sheet within the fundamental domain. We observe that $\mathcal{E}_{bend}$ has an inverse quadratic dependence on $R$, which signifies a Kirchhoff-Love type bending behavior(Reddy, 2006), i.e.,

[TABLE]

where $D$ can be interpreted as the modulus along the direction of bending, i.e., armchair or zigzag. This quadratic dependence on curvature is in good agreement with previous such studies for carbon nanotubes.Sánchez-Portal et al. (1999)

In Table 4, we list the bending modulus of the Xene sheets along the armchair and zigzag directions, obtained by fitting the data in Fig. 6. It is clear that the bending moduli of graphene are significantly larger than the bending moduli of the other Xenes. In particular, the bending modulus of graphene is factors of $3.3$ and $6.0$ larger than stanene in the zigzag and armchair directions, respectively. This is possibly a consequence of the short and strong bonds in graphene compared to the other Xenes, particularly stanene. We also find that, apart from graphene which is known to be close to isotropic Sánchez-Portal et al. (1999); Kudin et al. (2001); Wei et al. (2012), there is significant anisotropy in the bending modulus between the two directions. We correlate this anisotropy with the normalized buckled distance $\delta/a$, i.e., as $\delta/a$ increases, so does the anisotropy between the bending moduli along the two directions.

It is worth noting that the bending modulus of graphene computed here ($\sim 1.5$ eV) is in good agreement with previous such DFT predictions Kudin et al. (2001); Wei et al. (2012). The need for ab-initio calculations is clear from the scatter ($\sim 0.8-1.4$ eV Arroyo and Belytschko (2004); Lu et al. (2009)) in the predictions made while using empirical potentials . This need is further emphasized by the tremendously larger ($\sim 38.5$ eV Roman and Cranford (2014)) and therefore likely unphysical results obtained for the bending modulus of silicene using empirical potentials. Though state of the art efficient DFT implementations could possibly have been used to calculate the bending moduli of the Xenes studied in this work, the computational cost becomes prohibitively expensive, particularly as the radius of curvature approaches values representative of those found in experiments. Indeed, by exploiting the cyclic symmetry, the proposed formulation achieves a tremendous speedup, enabling the extremely efficient study of such systems, as quantified in the next subsection.

V.4 Scaling and performance

Finally, we turn to the scaling and performance of the proposed method. We choose zigzag silicon nanotubes as representative examples for this study and use discretization parameters $h=0.5$ Bohr and $N_{\eta}=3$, which results in energy and atomic forces that are converged to within $10^{-3}$ Ha/atom and $10^{-3}$ Ha/Bohr, respectively. These accuracies are more typical of those targeted in DFT calculations, including those involving geometry optimization and molecular dynamics. Indeed, as mentioned before, significantly more stringent accuracies were targeted in the previous two subsections to ensure that the scaling relations for the X nanotubes (i.e., bandgap as a function of the nanotube radius) and bending moduli of the Xene sheets were calculated to a high degree of precision.

We first perform a strong scaling study for a silicon nanotube with radius $6.1$ nm ($\mathfrak{N}=101$). Specifically, holding the system fixed, we increase the number of processors from $2$ to $152$ and determine the wall time associated with the complete simulation, i.e., total time for the calculation of the ground state electron density, energy, and atomic forces. We present the results so obtained in Fig. 7a, from which it is clear that we obtain good strong scaling, achieving an efficiency of $50\%$ on the largest number of processors relative to the smallest number of processors. The wall time on 152 processors is only 70 seconds, which is relatively small given the size of the system. Indeed, it is factor of $\mathcal{O}(100)$ smaller than SPARC, when run on the same number of cores 444For this system, since SPARC faces issues with SCF convergence when the mixing parameters used in this work are chosen, the reported speedup is estimated based on the time per SCF iteration.. SPARC itself has been shown to be significantly more efficient compared to established planewave codes like ABINIT, highlighting the efficiency of the current approach.

Next, we perform a weak scaling study by selecting a series of silicon nanotubes with radii from $6.1$ nm ($\mathfrak{N}=101$) to $100.06$ nm ($\mathfrak{N}=1650$), while correspondingly increasing the number of processors from $38$ to $619$. We choose $N^{proc}_{D}=1$ and set $N^{proc}_{K}$ such that each processor works on the symmetry-adapted Hamiltonians associated with $4$ discretized characters. We plot the SCF iteration time in Fig. 7b, from which it is clear that the proposed approach demonstrates good weak scaling for the range of systems and processors considered. Specifically, we obtain close to linear scaling, achieving $80\%$ efficiency for an increase in the nanotube radius by a factor of $\sim 16$. This can be justified by the fact that the work done per processor remains independent of system size in the proposed approach. Given that traditional DFT formulations instead scale cubically with the nanotube radius, systems such as the $100.06$ nm nanotube considered here would be exceedingly expensive, if not intractable, prior to this work.

It is worth noting that we have performed the weak scaling study by increasing the nanotube radius (i.e., $\mathfrak{N}$), while holding the system size within the fundamental domain fixed. Alternatively, we could have increased the system size within the fundamental domain, while holding the value of $\mathfrak{N}$ fixed. However, the results obtained in this case would be similar to those obtained by the underlying SPARC code Ghosh and Suryanarayana (2017a, b), and therefore are not reproduced here for the sake of brevity. It is also worth noting that, unlike in the strong scaling study where we report the wall time for the complete simulation, we report the wall time per SCF iteration for the weak scaling study. This is because of the increase in number of SCF iterations with the radius of the nanotube, a behavior that can be attributed to the change in electronic properties with system size.

The above results suggest that with sufficient computational resources, the proposed formulation allows for accurate DFT simulations of extremely large radius nanotubes, with modest wall times. To demonstrate this capability, we simulate a zigzag type III silicon tube of radius $R\sim 1$ $\mu$m ($\mathfrak{N}=16,473$) in $53$ minutes of wall time on just $353$ processors. We have found that the system has a negligible direct bandgap of $\mathcal{O}(10^{-4})$ eV at $\frac{\eta H}{2\pi}=0$, consistent with the results and scaling law obtained for silicon nanotubes in Section V.2, further verifying the accuracy of the simulation. Fig. 8 shows the band structure diagram for $\frac{\eta H}{2\pi}=0$, which bears a high degree of resemblance to the corresponding band structure diagram of the silicene sheet, as is to be expected, given the extremely large radius of the tube. To the best of our knowledge, this is the first example in literature of a fully resolved Kohn Sham DFT calculation of a system in which one of the length scales is of the order of a micron. Indeed, this example is slightly contrived since similar results could have been obtained by traditional DFT implementations at substantially reduced cost by using a flat sheet approximation of this nanotube, given its tremendously large radius and insignificant curvature induced effects. However, it serves to demonstrate the capabilities of the proposed method, with potential application to naturally large radius nanotubes Shin et al. (2004); McGary et al. (2006); Macak et al. (2008), where curvature induced effects are likely to be substantial.

In this work, though we have studied nanotube systems containing only 4 atoms in the fundamental domain, the proposed method is not restricted by this number, and given sufficient computational resources, the developed implementation can study systems containing up to a thousand atoms in the fundamental domain. Indeed, due to the significantly increased cost in such cases, the efficiency of the code would benefit from band parallelization as well as the AAR Pratapa et al. (2016); Suryanarayana et al. (2019) and DDBP Xu et al. (2018) methods. Even larger systems containing tens of thousands of atoms in the fundamental domain will then become accessible using the Complementary Subspace method Banerjee et al. (2018). Such advances will enable a number of applications, including the study of nanofilm bending Park et al. (2014); Haque and Saif (2003); Huang et al. (2005), which are intractable using traditional DFT methods.

VI Concluding remarks

In this work, we have developed a symmetry-adapted real-space formulation of Kohn-Sham DFT for cylindrical geometries and applied it to the study of large X (X=C, Si, Ge, Sn) nanotubes. Specifically, we have started from the original Kohn-Sham equations that are posed on all of space, and reduced them to the fundamental domain by accounting for the cyclic and periodic symmetries present in the angular and axial directions of the cylinder, respectively. We have implemented this approach for parallel computations using the high-order real-space finite-difference method, and verified its accuracy with respect to established planewave and real-space codes. We have used this implementation to study the band structure properties of X nanotubes and bending properties of Xene sheets. Specifically, we have first shown that zigzag and armchair X nanotubes with radii in the range of $1$ to $5$ nm are semiconducting, other than the armchair and zigzag type III carbon variants, for which we find a vanishingly small bandgap, indicative of metallic behavior. In particular, we have found that apart from armchair and zigzag type III carbon nanotubes, which demonstrate an inverse quadratic dependence of the bandgap with respect to radius, all other nanotubes demonstrate an inverse linear dependence. Next, we have exploited the the connection between cyclic symmetry and uniform bending deformations to calculate the bending moduli of Xene sheets in both zigzag and armchair directions for radii of curvature up to $5$ nm. We have found that the sheets obey Kirchhoff-Love type bending, with graphene and stanene demonstrating the largest and smallest moduli, respectively. In addition, apart from graphene, the sheets demonstrate significant bending anisotropy, with larger moduli along the armchair direction. Finally, we have shown that the proposed method is highly efficient and extremely well suited for parallel computations, which enables ab initio simulations of unprecedented size for systems with a relatively large degree of cyclic symmetry. In particular, we have shown that nanotubes with radii even at the micrometer scale can be simulated with modest computational resources and effort.

Overall, the proposed method provides an efficient framework for ab-initio simulations of 1D nanostructures with large radii as well as 1D/2D nanostructures under uniform bending, which are intractable using traditional formulations and implementations of DFT. This opens an avenue for the ab-initio study of the flexoelectric effect Hong and Vanderbilt (2013), in which a number of open questions remain Wang et al. (2019). The extension of the proposed method to include helical symmetry will enable the efficient ab-initio study of chiral nanotubes as well other nanostructures with helical symmetry, making it a worthy subject of future research.

Acknowledgements.

S.G. acknowledges support from the Army Research Laboratory which was accomplished under Cooperative Agreement Number W911NF-12-2-0022. A.S.B acknowledges support from the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, while at the Lawrence Berkeley National Laboratory. A.S.B also acknowledges support from the Minnesota Supercomputing Institute (MSI) for some of the computational resources that were used in this work. P.S. gratefully acknowledges the support of the National Science Foundation (CAREER-1553212). This research was supported in part through research cyberinfrastructure resources and services provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology, Atlanta, Georgia, USA. Some of the computations presented here were conducted on the Caltech High Performance Cluster partially supported by a grant from the Gordon and Betty Moore Foundation. The authors acknowledge the valuable comments and suggestions of the anonymous referees.