CARbon MAnifolds

Task D: Representing 3D Electron Densities on 2D Non-Euclidean Surfaces

Sub-Task D0: Developing a Finite Element / Finite Volumes code for Combinatorial Surfaces

In order to solve electronic structure equations on carbon manifolds (e.g. Kohn-Sham equations for electron density), the first step is to build a software library for solving differential equations on these non-Euclidean discrete surfaces. We will build a finite elements / finite volumes code that allow for numerical solution to partial differential equations directly on combinatorial surfaces (defined by the connectivity graph), without referring to an embedding in space. This will have applications in many areas of physics, as it provides a framework to numerically solve problems in discretizations of curved spaces, e.g. for numerical solutions in General Relativity, quantum gravity, or on the intricate gauge groups that appear in particle physics.

Two pilot MSc theses have laid the ground work for this task:

Simon K. Steensen, Wave Equations without Coordinates: Developing computational and mathematical methods for calculating wave equations without coordinates on two-dimensional manifolds arising from fullerene molecules. In it, Simon developed a first Finite Elements code for efficient solutions of differential equations on non-Euclidean manifolds, a fundamental component for solving wave equations for polyhedral molecules directly from topology/combnatoiral geometry without knowing how their 3D geometry will look. On top of this, a rudimentary 2D density functional theory was developed for interacting electrons moving along the manifold surface. However, this work was limited to linear elements on coarse meshes, limiting numerical accuracy.
Nikolai P. Nielsen's MSc, High-resolution Meshes for Two-dimensional Molecular Surfaces, adresses the problem of solving numerical differential equations accurately on non-Euclidean manifolds. The main challenge was to accurately discretize the Laplace-Beltrami operator under finite element mesh refinement, precomputing Gaussian curvature redistribution and locally remeshing the manifold to recover the refined metric. This work makes it possible to numerically accurate solve equations on carbon manifolds (and other non-Euclidean surface manifolds).

The lessons learned during this pilot project will form the starting point for the work by PhD1 and myself in Task D.

Sub-Task D1: Developing a Density Functional Theory for Fullerenes and Polyhedral Molecules

On top of the foundational work in D0, we will design a density functional theory (DFT) for non-Euclidean combinatorial surfaces. The main challenges are:

Electron densities cusp at atomic centres, which (if one does not design the representation carefully) requires extremely dense numerical representation near centres. This is simultaneously the least interesting and most difficult contribution to the electron density. We must design a representation in which the "plain graphene"-solution is factored out, and only the difference to this base is represented.
Develop exchange-correlation functionals for carbon manifolds. Since there is no global coordinate system, these should act from local information about electron densities (i.e., values and derivatives on individual or adjacent finite element cells), and local curvature of space. We will start with XC-functionals based on simple local density approximations (LDA), moving on to GGA methods, etc.
As we restrict the electron interaction to only act along the manifold surface (as we solve directly on the neighbour graph-induced metric): Can long-range inter-electron repulsion be meaningfully incorporated in the carbon manifold exchange-correlation functionals? Should the shortest geodesic distance be used, or contributions from all geodesic directions? Experiments comparing with full 3D DFT calculations with Gaussian (or other existing quantum chemistry software) can answer this question.

Sub-Task D2: Representing 3D Electron Densities in 2D and 1D

We do not just wish to capture an infinitely thin cross-section of the electron density along the surface. While we restrict electrons to move along the surface, the 3D distribution of electrons (normal to surface) in the real density is still important. How do we capture this information on the surface? That is: We want to represent a 3-dimensional constrained density on a 2-dimensional surface (on faces) and on 1-dimensional line segments (along bonds). One can think of this in analogy with holgraphic principles. We assume that electrons are centered around the surface and fall off as we move away (along the surface normal). Densities of bound electrons fall off exponentially with distance. Our ansatz shall therefore be:

Along edges (bonds) we assume that the density is approximately \begin{equation*} \rho_e(x,r) = a(x) e^{-b(x)r} \end{equation*} where $x$ is the coordinate along the edge, and $r = \sqrt{y^2+z^2}$ is the distance from the edge.
On the interior of the faces (pentagons and hexagons rings), we assume that the electron density is approximately \begin{equation*} \rho_f(x,y,r) = a(x,y) e^{-b(x,y)r} \end{equation*} where $x,y$ are 2D coordinates on the face-plane, and $r = |z|$ is the distance from the surface. As electron density is non-negative, $a(x,y) \ge 0$, and similarly $b(x,y) > 0$ must be strictly positive (or it would not describe a fall-off).

Edge Quadrature

Exponential fall-off away from edge: \begin{equation} \label{eq:edge-exp-integral} \begin{split} \int a(x) e^{-b(x)\sqrt{y^2+z^2}}~dx~dy~dz &= 2\pi \int_0^L dx \int_0^\infty dr~ a(x)\, r\, e^{-b(x)r}\\ &= 2\pi\int_0^L dx~\frac{a(x)}{b(x)^2}\\ &= 2\pi \sum_{q=1}^{n_q} \frac{a(x_q)}{b(x_q)^2} w_q \end{split} \end{equation} This means that we can do the 3D quadrature if we can make a 1D quadrature for rational functions. Alternatively, we can use a simple 1D Gauss-quadrature, but it is only approximate whenever $\frac{a(x)}{b(x)^2}$ is not a polynomial (we can use the least-squares polynomial fit).

Face Quadrature

Exponential fall-off away from face: $$\begin{eqnarray} \int_F a(x,y) e^{-b(x,y)|x|}~dx~dy~dz &= 2 \int_{(x,y)\in F} dx \int_0^\infty dr~ a(x,y)\, e^{-b(x,y) r}\\ &= 2 \int_{(x,y)\in F} dx~\frac{a(x,y)}{b(x,y)}\\ &= 2\sum_{q=1}^{N_q} \frac{a(x_q,y_q)}{b(x_q,y_q)} w_q \end{eqnarray}$$ This means that we can also do the 3D quadrature as a simple 2D rational-function quadrature for the faces!

Products of densities

$$\begin{eqnarray} \rho_1(x,y,z) \rho_2(x,y,z) &= a_1(x,y) e^{-b_1(x,y)|z|} a_2(x,y) e^{-b_2(x,y)|z|} \\ &= a_1(x,y) a_2(x,y) e^{-\left(b_1(x,y)+b_2(x,y)\right)|z|} \end{eqnarray}$$ I.e., the product is an exponential function with magnitude $a(x,y) = a_1(x,y) a_2(x,y)$ and fall-off exponent $b(x,y) = b_1(x,y) + b_2(x,y)$. Because of this, integrating face-densities over space can be evaluated as the 2D rational-function quadrature: $$\begin{eqnarray} \int_F \rho_1(x,y,z) \rho_2(x,y,z)~dx~dy~dz &= 2\sum_{q=1}^{N_q} w_q \frac{a_1(x_q,y_q) a_2(x_q,y_q)}{b_1(x_q,y_q)+b_2(x_q,y_q)} \end{eqnarray}$$ Since the $b$'s are strictly positive, the rational function has no roots in the integration domain, and efficient numerical quadratures can be computed.

Task:

This subtask will develop these ideas into efficient numerical codes that provide a finite-element formalism on polyhedral surfaces that faithfully represent the 3D electron densities investigated in Subtask D3 (described below) numerically in 2D. The ansatz of electron density falling off exponentially away from the molecular surface allows us to represent it on the surface with two 2D functions $a$ and $b$, and do 3D integration with 2D numerical rational-function quadratures.

Sub-Task D3: Studying Electron Densities of Fullerenes and Polyhedral Molecules

How do real electron densities (calculated by full ab initio quantum-chemistry with existing software) look for fullerenes? With optimized 3D geometries and full QC solutions for hundreds of fullerenes calculated at the highest feasible level of theory:

What are the actual electron density values sampled along the infinitesimally thin surface?
How does the real electron density fall off with distance orthogonal to the surface?
What is the best we can do with the 2D+exponential anzats representation of D2? This can be found with e.g. least-squares fits. Extrapolating this optimal 2D representation to space using the exponential ansatz of D2 back into 3D, how does the result compare to the full 3D density? What, if anything, important is missing? And can this be represented in 2D with a more elaborate ansatz?
Is my conjecture true, that the majority of the density will look simply like "a sheet of graphene wrapped around the polyhedral shape"? If not, can another simple systematic form be found, that we can factor out and treat without having to represent it on the numerical finite element mesh?
How well do the electron densities calculated with our newly developed DFT for non-Euclidean combinatorial geometry reproduce the real densities? Can we discover which physical effects we are missing that contribute importantly to the result? And can we discover a way to recover these within our computational framework?