# Point2CAD: Reverse Engineering CAD Models from 3D Point Clouds

Yujia Liu
ETH Zürich
Anton Obukhov
ETH Zürich
Jan Dirk Wegner
University of Zürich
Konrad Schindler
ETH Zürich

###### Abstract

Computer-Aided Design (CAD) model reconstruction from point clouds is an
important problem at the intersection of computer vision, graphics, and
machine learning; it saves the designer significant time when iterating on in-
the-wild objects. Recent advancements in this direction achieve relatively
reliable semantic segmentation but still struggle to produce an adequate
topology of the CAD model. In this work, we analyze the current state of the
art for that ill-posed task and identify shortcomings of existing methods. We
propose a hybrid analytic-neural reconstruction scheme that bridges the gap
between segmented point clouds and structured CAD models and can be readily
combined with different segmentation backbones. Moreover, to power the surface
fitting stage, we propose a novel implicit neural representation of freeform
surfaces, driving up the performance of our overall CAD reconstruction scheme.
We extensively evaluate our method on the popular ABC benchmark of CAD models
and set a new state-of-the-art for that dataset. Project page:
https://www.obukhov.ai/point2cad.

Figure 1:  Point2CAD is a pipeline for reconstructing complex CAD models,
including freeform surfaces, from 3D point clouds. A raw point cloud is
segmented into clusters corresponding to CAD faces, and each face is fitted
either with a geometric primitive or with a parametric freeform surface, using
a novel, implicit neural representation. Due to the analytic representation
the surfaces can be extended and intersected, such that topology emerges,
which is then used to clip the surface primitives.

##  1 Introduction

The task of reverse engineering CAD models from 3D point clouds has gained
increasing attention in recent years due to the rapid development of 3D
scanning technologies. Reverse engineering involves transforming a physical
object into a digital model in an editing-friendly format, which can be used
for analysis, visualization, and manufacturing.

Most approaches to reverse engineering CAD models follow a typical sequence of
steps. First, data capture and pre-processing are carried out to obtain point
cloud data sampled densely from the object’s surface. Subsequently, the point
cloud is segmented to identify individual regions of interest. The regions are
then classified to differentiate between different types of surfaces and
features. Analytical representations are generated based on the identified
regions, followed by the reconstruction of the model’s topology. Commercial
software systems offer interactive tools and functions for each step, such as
a dedicated segmentation or shape-fitting tool. Similarly, research that aims
to automate the process has focused chiefly on individual steps, whereas
little work covers the complete workflow of CAD model reconstruction.
Furthermore, most existing work is either limited to case studies of a few
models or is still quite far from producing satisfactory results.

In particular, reconstruction with simple parametric surface primitives has
been studied extensively. For such surfaces, both fitting and computing
intersections are comparatively straight-forward, and with sufficiently clean
point clouds, a boundary representation ($B$-rep) can be derived with
classical geometric methods [47, 1, 4]. On the contrary, while there are
several works about fitting freeform surfaces, often with various flavors of
splines, they mostly stop short of assembling individual surfaces into
complete $B$-reps [9, 42, 50]. In summary, geometrically complex free-form
surfaces and topologically complex assemblies of simple surfaces have mostly
been studied separately. The present work addresses this gap by proposing
Point2CAD, a method that recovers complete CAD models, including free-form
surfaces, edges, and corners. Refer to Sec. 7 of sup.mat. for notation.111A
note on terminology: many different terms have been used to denote the
topological elements of CAD models, see list in the supplementary material. We
use the terms _surface_ (2D), _edge_ (1D) and _corner_ (0D).

Our proposed pipeline for CAD model reconstruction from point clouds comprises
several steps. First, a pre-trained neural network segments the point cloud
into clusters corresponding to distinct surfaces. Second, surfaces are fitted
to the clusters, which involves both basic primitive fitting and a novel
neural fitting scheme for freeform surfaces. Third, adjacent surfaces are
intersected to recover edges, and adjacent edges are further intersected to
recover corners, thus obtaining a full $B$-rep. Taken together, these steps
form a comprehensive and versatile pipeline for reverse engineering point
clouds into CAD models; see Fig. 1. By combining modern, learning-based
segmentation backbones with classical geometric primitive fitting and with
recent neural field methods, we get the best of both worlds and obtain a
reconstruction pipeline that sets a new state of the art on the large and
diverse ABC benchmark [20]. Our contributions can be summarized as follows:

  * •

We propose Point2CAD, a comprehensive workflow for reconstructing complex
free-form parametric CAD models from raw 3D point clouds;

  * •

To power the surface fitting stage, we propose a novel neural representation
of freeform surfaces;

  * •

We demonstrate the superior quantitative and qualitative performance on the
challenging ABC dataset.

##  2 Related work

Our work aims to recover CAD models from raw point clouds. Towards this goal,
we draw inspiration from several domains of computer vision and 3D data
processing.

3D point cloud segmentation. Many studies segment point clouds by assigning
each point a class or instance label. Neural architectures to learn feature
embeddings of point clouds include: unordered sets of neighbors [33, 34],
graph convolutions [48, 23], point convolutions [25, 15], point-voxel learning
[37, 27], and transformers [52, 11]. Several authors use point-wise labeling
as a mechanism to structure point clouds into geometric primitives [12].
_E.g_., PrimitiveNet [16] segments point clouds into primitives via local
embedding and adversarial learning. In [28] graph convolution is used to
detect surface boundaries in 3D point clouds.

Primitive fitting [6, 18] is a common way to abstract 3D point clouds into
CAD-like parametric elements. Classical solutions rely on robust statistical
procedures like the Hough transform to recognize geometric primitives in point
clouds despite the noise and missing data [38]. Several learning-based
approaches have also been proposed to fit geometric primitives to point
clouds. SPFN [24] is trained to detect varying numbers of primitives at
different scales, supervised with ground truth surfaces and their membership
in the primitives. ParSeNet [42] finds parametric surfaces in point clouds,
including basic geometric primitives as well as B-spline surfaces, but does
not put much effort into connecting them. HPNet [50] concentrates on
partitioning the point cloud into segments using semantic and spectral
features and edge information in the form of an adjacency matrix but does not
fit actual primitives. Saporta et al. [40] proposed an unsupervised, recursive
neural architecture that learns to fit geometric primitives to 3D points,
while Vasu et al. [46] develop a hybrid representation that combines explicit
geometric shapes with signed distance functions to represent and manipulate
regular and freeform 3D shapes. Besides conventional geometric primitives,
parametric freeform surfaces are an important building block of CAD systems.
DeepSpline [9] utilizes a hierarchical RNN to reconstruct 2D spline curves, in
conjunction with an unsupervised learning approach to recover surfaces of
revolution and extrusion.

CAD representations. BRepNet [22] is a neural network architecture that
operates directly on $B$-rep data structures. UV-Net [17] introduces another
neural architecture that operates on $B$-rep data, without fitting to points.
JoinABLe [49] starts from CAD primitives and learns to assemble them and form
joints, using weak supervision available in standard CAD files. Another
popular CAD representation is constructive solid geometry (CSG). CAPRI-Net
[51] learns to fit a collection of quadric 3D primitives, and an associated
CSG-tree, while CSGNet [41] learns to generate CSG programs for given 2D or 3D
shapes, using a convolutional encoder and a recurrent decoder.

Generic CAD modeling. Reconstructing a complete CAD model [31, 3] is
challenging, as the output should include both the geometric surface shapes
and their connectivity. Traditional methods for reverse engineering often
employ iterative RANSAC estimation to extract primitive shapes from point
clouds, _e.g_., [26]. Benko et al. [2] outline a hand-engineered, sequential
procedure to reconstruct CAD models via segmentation, surface and topology
fitting, and blending. Recent, learning-based approaches in the same spirit
include Point2Cyl [45], a neural network that transforms a 3D point cloud into
a set of extrusion cylinders by predicting point segmentation, base labels,
normals, and extrusion parameters in closed form. ComplexGen [10] reconstructs
CAD models by autoregressively detecting geometric primitives. The model
consists of a sparse CNN encoder, three transformer decoders for geometric
primitives and topology, and post-processing with global optimization to
refine the rather rough transformer predictions.

Manifold learning refers to a family of methods that aim to discover a low-
dimensional manifold underlying a higher-dimensional, possibly noisy, data
set. Manifold learning methods can be understood as non-linear extensions of
PCA [7] attempting to “flatten” curved manifolds. Most methods focus on
visualization or fidelity to particular input points, by explicitly using
those points to construct the mapping [39, 44, 29]. Autoencoders provide a
natural way to not only project data non-linearly to a lower-dimensional
latent space but also to decode back from the latent space to the original
data space [21]. We extend this approach with the recent findings about
implicit neural representations to design a fitting method for freeform
surfaces.

Implicit Neural Representations (INR) are a framework to encode an arbitrary
function, often observed in the form of sparse samples, into the weights of a
neural network. Recent research about neural rendering has led to several
useful and transferrable insights about implicit neural fields, for instance,
the importance of positional encoding [30] and new (e.g., sinusoidal)
activation functions [43]. These findings form the basis for our freeform
surface fitting method.

Figure 2:  Visualization of open (left) and closed (right) surface fitting
with the proposed INR. Top: point cloud sampled from a freeform surface;
middle: INR fit, showing high or low principal curvature; bottom: same fit
with color indicating $uv$ coordinates. Our INR reproduces curvature where
necessary and smoothly extrapolates without oscillations to facilitate
subsequent steps.

Figure 3:  Inspiration items of our INR surface fitting. UMAP [29] (1st
column, middle) learns the underlying 2D manifold of 3D points along with the
inverse mapping, but cannot capture its smooth analytic representation. Early
experiments with ReLU activations confirmed its low-frequency bias (1st
column, bottom) and piece-wise linearity. SiLU [13] (2nd column) suffers from
low-frequency bias, too. We found that resolving it by adding positional
encoding is challenging under our training protocol (3rd column). SIREN [43]
(4th column) fits the data precisely, however, does not extrapolate well. A
combination of the best parts of these ingredients led us to our INR surface
fitting method, seen in Fig. 2 and discussed in Sec. 3.2.

##  3 Method

Contrary to the recent trend towards generic, end-to-end deep learning
pipelines, we found it advantageous to split the reconstruction process into
steps and only use neural methods where necessary. Overall, our method
consists of the following stages, _cf_. Fig. 1:

  1. 1.

Partition the point cloud into clusters corresponding to the CAD model’s
topological faces. We rely on existing (pretrained) neural network methods for
that step.

  2. 2.

Fit an analytical surface primitive to each cluster. Here we use a hybrid
approach: First, we test a set of prevalent geometric surface primitives that
admit efficient closed-form fitting. We propose a novel fitting scheme for
freeform surfaces based on neural representations.

  3. 3.

Identify the effective area of each (finite) parametric surface and clip it,
leaving enough margin to ensure intersections of adjacent surfaces.

  4. 4.

Perform pairwise surface intersections to obtain a set of topologically
plausible object edges. Using these edges, remove parts not supported by input
points.

  5. 5.

Perform pairwise edge intersection to identify a set of topological corners.
Clip edges based on proximity to the remaining surface regions using these
corners.

As a result of applying Point2CAD, we obtain a CAD model in $B$-rep format,
which includes analytical surfaces to represent the model’s faces, compatible
edges and corners, and adjacency matrices that encode topology.

A cornerstone of the proposed pipeline is a novel, efficient method for
fitting the implicit neural representation of a freeform, parametric
2D-manifold surface in 3D to an unordered set of points. A particular
difficulty in this context is that the surface must be extrapolated beyond the
supporting points to enable robust surface intersections, an ability that
tends to be challenging for neural surface estimators. To that end, we design
a novel freeform surface parameterization, drawing inspiration from recent
advances in neural fields, positional encoding schemes, and efficient test-
time optimization. See Sec. 3.2 for details.

While using neural methods, INR is a pure test-time optimization. We
effectively side-step learning priors from training data except for the
initial low-level segmentation. We argue that simple analytical procedures
like surface fitting do not necessarily benefit from learned priors and that
an over-reliance on data-driven learning may have hampered recent work on
reverse engineering CAD models. Moreover, separately predicting – rather than
constructing – the topology does not guarantee that it is consistent with the
geometry, _e.g_., edges could be instantiated such that they do not lie in a
surface, see Fig. 5. We empirically support this claim in our experiments,
where we outperform purely learning-based methods like [10], which are
inherently plagued by even slight differences between the training and test
distributions. When designing Point2CAD, we aimed for a generalizable scheme
that uses offline learning where necessary – in particular, for feature
learning and segmentation, ill-posed tasks that are notoriously hard to
engineer heuristically; but that does not replace well-understood analytical
operations with less robust, data-driven feed-forward modeling. As a side
effect, obviating parameter-heavy neural networks, where possible, also
results in a computationally less demanding method.

###  3.1 Parametrization of standard primitives

The geometric primitives that our current implementation handles are planes,
spheres, cylinders, and cones. The exact parameterizations are taken as in
[42], see Sec. 8 of sup.mat.

We observed that utilizing the predicted surface type may harm the
reconstruction since the prediction is not always correct. Therefore, we
circumvent that step and instead rely on an exhaustive, data-driven
comparison: we fit each primitive type (including INR) to point clusters and
select the simplest model with the lowest reconstruction error.

###  3.2 Freeform surface parametrization with INR

In the context of the overall CAD reconstruction pipeline, this step should
satisfy several requirements: (1) resilience to noise in the input point
cloud, which is inevitable with real 3D sensors; (2) support for inverse
mapping, to enable traversal of the latent space; (3) flexibility in
interpolation mode, to ensure high data fidelity and avoid over-smoothing; (4)
strong regularization in extrapolation mode, required to construct smooth,
plausible margins around the observed data points for the subsequent surface
intersection; (5) fast fitting with low computational cost.

Existing manifold learning techniques that support the inverse transform, such
as UMAP [29], could not be readily used for this task, as demonstrated in Fig.
3 (1st column, middle). We thus developed a custom neural autoencoder [21]
with a single hidden layer, which is sufficient for simple non-linear
transforms [14].

We found that for our purposes, the choice of activation function is critical:
Besides causing a non-smooth, piece-wise linear approximation, ReLU
activations [8] exhibit a noticeable low-frequency bias [35], as can be seen
in Fig. 3 (1st column, bottom). High-frequency undulations are over-smoothed,
such that entire point cloud regions noticeably deviate from the fit on either
side. Replacing ReLU with SiLU [13] enhances surface smoothness, but does not
overcome the low-frequency bias (2nd column). Adding positional encoding [30]
in conjunction with SiLU resolves the issue, as one would expect from the
literature, but provokes new artifacts, due to the fixed set of spatial
frequencies (3rd column). Sinusoidal activations [43] (4th column) achieve
excellent data fitting, but impose little regularization when not supported by
data, and thus extrapolate poorly.

A combination of activations with desired properties within the same layer
overcomes these individual problems: equipping some neurons with SiLU and the
remaining ones with sinusoidal activations produces geometrically faithful
fits that extrapolate smoothly, see Fig. 2. Remarkably, this behavior emerges
without any dedicated mechanisms to steer the collaboration between the two
types of neurons.

Figure 4:  Block scheme of our proposed INR for freeform surface fitting. We
encode groups of 3D surface points into a latent 2D $uv$ space. We use a
mixture of activations to achieve high fidelity of fitting and smooth
extrapolation. Both open and closed surface fitting is supported via
preconfigured routing of activations in the latent space. See Sec. 3.2 for
more details.

Fig. 4 depicts a block diagram of the proposed autoencoder. Each surface is
autoencoded independently, by feeding batches of its 3D points ($X$,$Y$,$Z$)
through a 1-layer MLP encoder and a corresponding decoder, with a 2D
bottleneck corresponding to manifold coordinates ($u$,$v$). The 2D latent
space is regularized to a unit square $[-1,1]\times[-1,1]$. Both closed, and
open surfaces (independently in the $u$ and $v$ dimensions) are supported via
pre-configured routing.

For each surface, the weights of a template INR are initialized randomly and
optimized with standard mini-batch (or full-batch) descent. We run Adam [19]
for 1000 steps, with 50 steps warm-up of the learning rate, followed by a
linear decay that reaches zero at the last step. The complete optimization
takes only a few seconds on a single GPU, and multiple surfaces can be fitted
in parallel.

For extrapolation and latent space traversal, we encode all cluster points
into the latent $uv$ space and store the bounding box parameters along with
the autoencoder. We reset the corresponding axis range for surfaces with
closed dimensions to $[-1,1]$. To sample the extended surface with the margin,
we extend the bounding box by $10\%$ in both dimensions and compute 3D surface
points using the decoder.

###  3.3 Topology reconstruction

Following the instantiation of all individual surfaces, it becomes necessary
to establish their boundaries, which amounts to finding the intersection
curves between adjacent surfaces. Since this step is crucial for subsequent
stages of CAD reconstruction, it is important to correctly identify and
recover the intersection curves. This requirement again calls for
deterministic geometric methods, as independent detection can hardly guarantee
a set of complete edges consistent with the surfaces.

While it is theoretically possible to work out the surface intersections
analytically, multiple freeform surfaces lead to extremely complicated, and
potentially unstable calculations. Hence, we opt for a generic solver and
convert all surfaces to triangle meshes for this step to be accurately
intersected with mature, numerically stable computational tools. For each
surface, the (infinite) geometric primitive is trimmed to form a margin of
width $\epsilon$ around the supporting points to ensure enough space for the
subsequent steps and then fed to a standard meshing algorithm by tessellation
and then triangulation. We identify all intersecting surfaces for each meshed
surface, compute the pairwise mesh intersections to obtain poly-line edges,
and re-mesh the surfaces along those edges. Similarly, we intersect adjacent
edge poly-lines to obtain corner points and, in turn, use those corners to
trim the edges (which still extend beyond the surfaces since they were
generated from unclipped surfaces). See Sec. 9 of sup.mat. for a detailed
algorithm.

Such a discretization step temporarily side-steps the $B$-rep format. However,
it is possible to compute the true analytical intersections starting from the
discrete approximations. This arguably simpler task is left for future work.

Figure 5:  The final results gallery, with different colors denoting different
topological surfaces, and edges and corners depicted in black. The first
three, left to right: input point cloud, ground truth mesh, and reconstruction
with ComplexGen [10]. The last three are the proposed Point2CAD method applied
to different segmentation of the input points: with HPNet [50], ParSeNet [42],
and Ground Truth. Our method reconstructs the ground truth geometry and
topology from the ground truth segmentation nearly perfectly. When applied on
top of pretrained segmentation modules, it outperforms the competition by a
high margin.

##  4 Experiments

The rather flexible constraints on the input to our method permit employing it
in various setups. Such may include:

  1. –

Evaluation of the ground truth point cloud clustering or segmentation as a way
to quantify the contribution of point cloud sampling sparsity and sample noise
on the reconstruction quality (aliased “Point2CAD GT”);

  2. –

Usage of geometry fitting and topology extraction on top of any pretrained
point cloud clustering or segmentation methods, such as ParSeNet [42] or HPNet
[50] (aliased “Point2CAD $\langle$method$\rangle$”).

Under both setups, we first individually evaluate the surface reconstruction
quality of the most challenging freeform surfaces, which is one of the main
contributions of our work. We then analyze the geometry and topology of the
entire reconstructed CAD dataset. This enables us to provide an overall
assessment of the effectiveness of our approach in reconstructing surfaces and
CAD models from point clouds.

###  4.1 Dataset

To evaluate the performance of Point2CAD, we conduct experiments on the ABC
dataset [20], a large-scale collection of CAD models ($\sim$1.000.000 models)
widely used in geometric deep learning. We use the same split as ParseNet,
where each model contains at least one freeform surface, to facilitate a fair
comparison with existing methods and to better demonstrate the feasibility of
our approach.

###  4.2 Evaluation metrics

  * –

Residual error measures the discrepancy between the reconstructed surface and
the corresponding ground truth surface obtained by Hungarian matching.
Following the established protocol [42, 50], we densely sample points on the
ground truth surface and evaluate their correspondence with the reconstructed
surface. Specifically, we compute the distance between each sampled point and
its projected point on the reconstructed surface, which enables us to assess
how well the fitting algorithm has captured the underlying geometry of the
original surface. The residual error for $i$-th surface is:

$Err_{i}=\frac{1}{N_{i}}\sum_{\mathbf{s}_{i,k}\in
S_{i}^{\text{gt}}}d(\mathbf{s}_{i,k},S_{i}^{r}),$ (1)

where $N_{i}$ is the number of sampled points, $\mathbf{s}_{i,k}$ is the
$k$-th sampled point on the $i$-th ground truth surface $S_{i}^{\text{gt}}$,
and its corresponding reconstructed surface is $S_{i}^{r}$. The residual error
$Err_{i}$ is computed as the average Euclidean distance $d(\cdot,\cdot)$
between the sampled points and their corresponding projections on the
reconstructed surface. The residual error for $i$-th surface is a simple
average can provide an overall assessment of a CAD model,
$Err=\frac{1}{M}\sum_{i=1}^{M}Err_{i}$.

  * –

P-coverage. It quantifies the proportion of input point cloud covered by the
generated surface:

$Pcov=\frac{1}{|P|}\sum_{k=1}^{|P|}\mathbb{I}\\{d(p_{k},S^{r})\leq r\\},$ (2)

where $|P|$ is the number of points in the entire input point set $P$ and $S$
is all CAD reconstruction surfaces. $r$ is a distance threshold, set as $0.01$
following previous works.

  * –

Surface precision, recall, F-score. We adopted several evaluation metrics
introduced in ComplexGen [10] to assess the efficacy of our method in
generating reconstructed surfaces. A surface is considered a true positive if
the distance between it and its closest ground truth surface is below a
certain threshold $\theta$, which serves as a measure of accuracy for the
method in a geometric sense.

  * –

Edge precision, recall, F-score. The effectiveness of the proposed method can
also be assessed from a topological perspective, by computing the precision,
recall, and F-score on the edges resulting from pairwise surface
intersections. This evaluation strategy provides a reliable measure of the
method’s ability to capture the topology of the reconstructed surfaces.

  * –

Corner precision, recall, F-score. Likewise, an accurate estimation of corner
positions resulting from edge intersections is an essential indicator of
overall accuracy.

  * –

Chamfer distance. This metric assesses the accuracy of the reconstructed CAD
model as a whole, without distinguishing between individual surfaces. To that
end, we sample the entire surface of both the reconstructed model and the
ground truth object and calculate the chamfer distance between the two point
sets.

###  4.3 Evaluation on freeform surfaces

Both ComplexGen [10] and ParseNet [42] have successfully demonstrated the
prediction of the freeform surfaces. Both papers employ a custom
parameterization that represents any freeform surface as a structured grid of
20x20 control points, regardless of the underlying scale.

Table 1:  Comparison of surface fitting. ComplexGen uses a transformer decoder
to convert a surface latent code to a 20x20 control-points grid. ParseNet
employs a neural network for the same kind of output. Our method utilizes the
proposed INR surface format.

| Open surfaces | Closed surfaces
---|---|---
| Res-err $\downarrow$ | P-cover $\uparrow$ | Res-err $\downarrow$ | P-cover $\uparrow$
ComplexGen | $0.021$ | $0.938$ | $0.023$ | $0.900$
ParseNet | $0.006$ | $0.930$ | $0.008$ | $0.902$
Ours | $\mathbf{0.002}$ | $\mathbf{0.999}$ | $\mathbf{0.003}$ | $\mathbf{0.998}$

Table 2: Ablation of INR design performance with and w/o noise added to the
input points. Our method demonstrates favorable properties as measured by
residual error and P-coverage.

|  | Open surfaces | Closed surfaces
---|---|---|---
Noise | Activation | Res-err $\downarrow$ | P-cover $\uparrow$ | Res-err $\downarrow$ | P-cover $\uparrow$
$\sigma=0.05$ | ReLU | $0.016$ | $0.519$ | $0.018$ | $0.509$
SiLU | $0.024$ | $0.522$ | $0.030$ | $0.543$
SiLU+PosEnc | $0.014$ | $0.561$ | $0.015$ | $0.516$
SIREN | $0.012$ | $0.524$ | $0.014$ | $0.521$
Ours | $\mathbf{0.011}$ | $\mathbf{0.585}$ | $\mathbf{0.013}$ | $\mathbf{0.581}$
No noise | ReLU | $0.006$ | $0.986$ | $0.008$ | $0.981$
SiLU | $0.020$ | $0.905$ | $0.023$ | $0.887$
SiLU+PosEnc | $0.004$ | $0.992$ | $0.007$ | $0.982$
SIREN | $0.003$ | $0.996$ | $0.006$ | $0.990$
Ours | $\mathbf{0.002}$ | $\mathbf{0.999}$ | $\mathbf{0.003}$ | $\mathbf{0.998}$

Table 3: Ablation of INR under surface margin extrapolation by a fraction $c$
of model diameter. Point2CAD uses $c=10\%$.

Res-err $\downarrow$ | $c=0\%$ | $c=10\%$ | $c=20\%$
---|---|---|---
SiLU | $0.0582$ | $0.0565$ | $0.0714$
SIREN | $0.0195$ | $0.0376$ | $0.1205$
Ours | $\mathbf{0.0052}$ | $\mathbf{0.0156}$ | $\mathbf{0.0474}$

Table 4: Geometric evaluation of reconstructed CADs. Segmentation denotes the
method used for point cloud clustering. “GT” stands for oracle ground truth
segmentation, which is also an upper bound of the performance of our method.

| Segmentation | Res-err $\downarrow$ | P-cover $\uparrow$ | Chamfer $\downarrow$
---|---|---|---|---
ComplexGen | N/A | $0.020$ | $\mathbf{0.950}$ | $0.042$
Point2CAD | ParseNet | $\mathbf{0.018}$ | $0.942$ | $\mathbf{0.017}$
Point2CAD | HPNet | $0.020$ | $0.937$ | $0.018$
Point2CAD | GT | $0.011$ | $0.968$ | $0.016$

To allow for an apples-to-apples comparison, we seek to obtain the results of
all three methods on distinct point cloud clusters. While this is
straightforward with ParseNet and our method, ComplexGen does not operate on
individual clusters but instead requires the entire point cloud at once and
outputs all $B$-rep elements at once, too. Thus, to obtain ComplexGen
predictions for a particular cluster, we run it with the entire point cloud
input, find the nearest predicted surface to the queried cluster out of all
the predictions, and perform the surface-to-surface evaluation on that pair.

We repeat this procedure for a subset of freeform surfaces of the ABC dataset
and average the geometry metrics; see Tab. 1 for quantitative evaluation. As
can be seen, our INR fitting method generates freeform surfaces that represent
the underlying points more faithfully than prior parameterizations.
Qualitative results of Point2CAD freeform surface fitting confirm the
qualitative study, which can be seen in Fig. 2, and the side-by-side
comparison is in Fig. 5.

Additionally, we studied the behavior of different INR activations under
noise-induced corruption of the input point clouds (Tab. 2) and surface margin
extrapolation (Tab. 3), a crucial property of our pipeline. For the latter,
given a point cloud patch of a freeform surface, we select a subset of points
from the center of the patch for fitting, leaving a margin of $c\%$ for
extrapolation evaluation. Our INR design consistently performs well across
various settings.

###  4.4 Reconstructed CAD evaluation

We conduct an evaluation of Point2CAD in terms of the reconstruction quality
for surfaces, edges, and corners. We compare our results to those obtained
with ComplexGen [10], which is also capable of generating all those geometric
elements and their topological relationships. However, it should be noted that
the models generated by ComplexGen are often invalid, in the sense that the
predicted geometry and topology are inconsistent. We are not aware of an easy
way to rectify this issue. A possible solution could be only to utilize the
surfaces generated by ComplexGen and apply our method to refit and intersect
them.

Table 5: Evaluation on CAD Reconstruction of Surfaces, Edges and Corners from
the aspects of accuracy and completeness. The surface evaluation assesses the
reconstruction performance in geometric terms, while these metrics on edges
and corners reflect the performance of reconstruction in terms of topological
properties. Given that the distance calculation is based on sampling, we adopt
various threshold values for evaluating surfaces, edges and corners.
Specifically, the threshold for a single-point corner is expected to be
smaller than that of edges, while the threshold for edges should be smaller
than that of surfaces, which aligns with the intuitive expectations in CAD
evaluation.

Surfaces | $\theta_{\text{surface}}=0.08$ | $\theta_{\text{surface}}=0.06$ | $\theta_{\text{surface}}=0.03$
---|---|---|---
Method | Segmentation | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$ | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$ | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$
ComplexGen | N/A | $0.732$ | $0.732$ | $0.731$ | $0.633$ | $0.641$ | $0.637$ | $0.370$ | $0.388$ | $0.379$
Point2CAD | ParseNet | $0.817$ | $\mathbf{0.750}$ | $0.782$ | $0.764$ | $0.693$ | $0.727$ | $0.624$ | $0.562$ | $0.591$
Point2CAD | HPNet | $\mathbf{0.836}$ | $0.745$ | $\mathbf{0.788}$ | $\mathbf{0.796}$ | $\mathbf{0.695}$ | $\mathbf{0.742}$ | $\mathbf{0.697}$ | $\mathbf{0.593}$ | $\mathbf{0.641}$
Point2CAD | GT | $0.976$ | $0.920$ | $0.947$ | $0.962$ | $0.898$ | $0.929$ | $0.873$ | $0.801$ | $0.836$
Edges | $\theta_{\text{edge}}=0.05$ | $\theta_{\text{edge}}=0.03$ | $\theta_{\text{edge}}=0.02$
Method | Segmentation | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$ | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$ | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$
ComplexGen | N/A | $0.620$ | $0.576$ | $0.597$ | $0.421$ | $0.397$ | $0.409$ | $0.291$ | $0.279$ | $0.285$
Point2CAD | ParseNet | $0.636$ | $\mathbf{0.596}$ | $\mathbf{0.615}$ | $0.523$ | $\mathbf{0.486}$ | $\mathbf{0.504}$ | $0.440$ | $\mathbf{0.414}$ | $\mathbf{0.427}$
Point2CAD | HPNet | $\mathbf{0.637}$ | $0.586$ | $0.611$ | $\mathbf{0.526}$ | $0.473$ | $0.498$ | $\mathbf{0.449}$ | $0.397$ | $0.421$
Point2CAD | GT | $0.863$ | $0.774$ | $0.816$ | $0.766$ | $0.673$ | $0.717$ | $0.686$ | $0.597$ | $0.639$
Corners | $\theta_{\text{corner}}=0.03$ | $\theta_{\text{corner}}=0.02$ | $\theta_{\text{corner}}=0.01$
Method | Segmentation | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$ | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$ | precision $\uparrow$ | recall $\uparrow$ | F-score $\uparrow$
ComplexGen | N/A | $\mathbf{0.667}$ | $0.633$ | $0.650$ | $0.483$ | $0.454$ | $0.468$ | $0.217$ | $0.203$ | $0.210$
Point2CAD | ParseNet | $0.661$ | $\mathbf{0.641}$ | $\mathbf{0.651}$ | $\mathbf{0.553}$ | $\mathbf{0.529}$ | $\mathbf{0.541}$ | $\mathbf{0.398}$ | $\mathbf{0.392}$ | $\mathbf{0.395}$
Point2CAD | HPNet | $0.646$ | $0.613$ | $0.629$ | $0.549$ | $0.521$ | $0.535$ | $0.389$ | $0.386$ | $0.388$
Point2CAD | GT | $0.780$ | $0.696$ | $0.736$ | $0.704$ | $0.662$ | $0.661$ | $0.581$ | $0.515$ | $0.546$

As seen in Tab. 4, Point2CAD outperforms ComplexGen in most metrics.
Furthermore, Point2CAD achieves excellent results when fed with ground truth
segmentations, indicating that it will likely further improve as better point
cloud segmentation engines become available. We also observed that our
P-coverage results were slightly inferior to ComplexGen’s. This can be
attributed to the fact that, in the final step of our reconstruction, we trim
the margins of the individual surfaces to establish a clear and well-defined
boundary. As a result, the coverage of our reconstructed CAD models is
somewhat reduced compared to ComplexGen, which does not perform this step and
typically generates more surfaces than the ground truth has. The presence of
numerous spurious surfaces naturally inflates the P-coverage. This limitation
of the metric is confirmed by the fact that we obtain higher P-coverage if we
evaluate Point2CAD without clipping the surfaces at the edges
(Point2CAD+ParseNet: $0.956$, Point2CAD+HPNet: $0.953$, Point2CAD+GT:
$0.971$).

We also employ several standard evaluation metrics, including precision,
recall, and F-score, to measure the accuracy and completeness of the surfaces
after the overall reconstruction of CAD under different distance thresholds
$\theta_{\text{surface}}$, see Tab. 5 (top). The results suggest that
Point2CAD with HPNet segmentation backbone achieves the highest overall
performance, while reconstruction based on the ground truth point cloud
segmentation provides the best results, as expected.

We utilize the precision, recall, and F-score of the reconstructed curves and
corners as metrics to evaluate the topological reconstruction performance;
results are shown in Tab. 5 (middle, bottom).

##  5 Conclusion

We present Point2CAD, a comprehensive and versatile approach for reverse
engineering CAD models from point clouds, which can yield various types of
surfaces and topologically correct reconstructions. The proposed method
segments point clouds with a pre-trained segmentation backbone but then
employs learning-free optimization methods to fit geometric primitives,
including freeform surfaces that are optimized with a novel neural surface
estimation scheme. The recovered surfaces are analytically intersected to
obtain the edges and corners of the model. Empirically, our proposed method
qualitatively and quantitatively outperforms prior art on the popular ABC
dataset of CAD models. In our view, these results support our hypothesis that
the main role of machine learning in the context of CAD reconstruction is
actually segmentation, whereas subsequent steps are adequately solved by
classical geometric operations and may not always require learning.

Limitations While our method offers a solution to CAD model reconstruction, it
is subject to certain limitations. Notably, our method is sensitive to the
quality of the point cloud segmentation. A number of outliers in the
prediction of points belonging to the same surface can adversely impact the
subsequent surface fitting and topology reconstruction. Despite the fact that
our method is not end-to-end-trainable, the main guideline in our design is to
divide the pipeline into isolated, accountable steps and only employ neural
methods where needed. While this runs against a puristic deep learning spirit,
our modular pipeline outperforms the end-to-end learned solutions. We see this
as an important message, although we do not rule out that it could serve as a
foundation for future end-to-end works.

Future work One promising direction would be to exploit the proposed INR
fitting scheme in a feedback loop of any other algorithm to drive point cloud
partitioning into clusters, thus lifting the dependence of the reconstruction
pipeline on the quality of point cloud segmentation. Differently, bringing
back learned priors and the end-to-end property should make it possible to
further improve the results.

## References

  * Alai [2013] Shashank Alai.  A review of 3d design parameterization using reverse engineering.  _International Journal of Emerging Technology and Advanced Engineering_ , 3(10):171–179, 2013.
  * Benkő et al. [2001] Pál Benkő, Ralph R Martin, and Tamás Várady.  Algorithms for reverse engineering boundary representation models.  _Computer-Aided Design_ , 33(11):839–851, 2001\.
  * Buonamici et al. [2018] Francesco Buonamici, Monica Carfagni, Rocco Furferi, Lapo Governi, Alessandro Lapini, and Yary Volpe.  Reverse engineering modeling methods and tools: a survey.  _Computer-Aided Design and Applications_ , 15(3):443–464, 2018.
  * Durupt et al. [2008] Alexandre Durupt, Sebastien Remy, Guillaume Ducellier, and Benoît Eynard.  From a 3d point cloud to an engineering CAD model: a knowledge-product-based approach for reverse engineering.  _Virtual and Physical Prototyping_ , 3(2):51–59, 2008.
  * Eberly [2006] David Eberly.  Geometric tools, 2006.
  * Farin [2002] Gerald Farin.  _Curves and surfaces for CAGD: a practical guide_.  Morgan Kaufmann, 2002.
  * F.R.S. [1901] Karl Pearson F.R.S.  Liii. on lines and planes of closest fit to systems of points in space.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 2(11):559–572, 1901.
  * Fukushima [1975] Kunihiko Fukushima.  Cognitron: A self-organizing multilayered neural network.  _Biological cybernetics_ , 20(3-4):121–136, 1975\.
  * Gao et al. [2019] Jun Gao, Chengcheng Tang, Vignesh Ganapathi-Subramanian, Jiahui Huang, Hao Su, and Leonidas J Guibas.  DeepSpline: Data-driven reconstruction of parametric curves and surfaces.  _arXiv preprint arXiv:1901.03781_ , 2019.
  * Guo et al. [2022] Haoxiang Guo, Shilin Liu, Hao Pan, Yang Liu, Xin Tong, and Baining Guo.  ComplexGen: CAD reconstruction by B-rep chain complex generation.  _ACM TOG_ , 41(4):1–18, 2022.
  * Guo et al. [2021] Meng-Hao Guo, Jun-Xiong Cai, Zheng-Ning Liu, Tai-Jiang Mu, Ralph R Martin, and Shi-Min Hu.  PCT: Point cloud transformer.  _Computational Visual Media_ , 7:187–199, 2021.
  * He and Wang [2018] Chen He and Chunmeng Wang.  A survey on segmentation of 3d models.  _Wireless Personal Communications_ , 102:3835–3842, 2018\.
  * Hendrycks and Gimpel [2016] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus).  _arXiv preprint arXiv:1606.08415_ , 2016.
  * Hornik et al. [1989] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  _Neural networks_ , 2(5):359–366, 1989.
  * Hua et al. [2018] Binh-Son Hua, Minh-Khoi Tran, and Sai-Kit Yeung.  Pointwise convolutional neural networks.  In _CVPR_ , 2018.
  * Huang et al. [2021] Jingwei Huang, Yanfeng Zhang, and Mingwei Sun.  PrimitiveNet: Primitive instance segmentation with local primitive embedding under adversarial metric.  In _ICCV_ , 2021.
  * Jayaraman et al. [2021] Pradeep Kumar Jayaraman, Aditya Sanghi, Joseph G Lambourne, Karl DD Willis, Thomas Davies, Hooman Shayani, and Nigel Morris.  UV-Net: Learning from boundary representations.  In _CVPR_ , 2021.
  * Kaiser et al. [2019] Adrien Kaiser, Jose Alonso Ybanez Zepeda, and Tamy Boubekeur.  A survey of simple geometric primitives detection methods for captured 3d data, 2019.
  * Kingma and Ba [2015] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.
  * Koch et al. [2019] Sebastian Koch, Albert Matveev, Zhongshi Jiang, Francis Williams, Alexey Artemov, Evgeny Burnaev, Marc Alexa, Denis Zorin, and Daniele Panozzo.  ABC: A big CAD model dataset for geometric deep learning.  In _CVPR_ , 2019.
  * Kramer [1991] Mark A Kramer.  Nonlinear principal component analysis using autoassociative neural networks.  _AIChE journal_ , 37(2):233–243, 1991.
  * Lambourne et al. [2021] Joseph G Lambourne, Karl DD Willis, Pradeep Kumar Jayaraman, Aditya Sanghi, Peter Meltzer, and Hooman Shayani.  BRepNet: A topological message passing system for solid models.  In _CVPR_ , 2021.
  * Li et al. [2019a] Guohao Li, Matthias Muller, Ali Thabet, and Bernard Ghanem.  DeepGCNs: Can GCNs go as deep as CNNs?  In _ICCV_ , 2019a.
  * Li et al. [2019b] Lingxiao Li, Minhyuk Sung, Anastasia Dubrovina, Li Yi, and Leonidas J Guibas.  Supervised fitting of geometric primitives to 3d point clouds.  In _CVPR_ , 2019b.
  * Li et al. [2018] Yangyan Li, Rui Bu, Mingchao Sun, Wei Wu, Xinhan Di, and Baoquan Chen.  PointCNN: Convolution on x-transformed points.  _NeurIPS_ , 2018.
  * Liu [2020] Jin Liu.  An adaptive process of reverse engineering from point clouds to CAD models.  _International Journal of Computer Integrated Manufacturing_ , 33(9):840–858, 2020.
  * Liu et al. [2019] Zhijian Liu, Haotian Tang, Yujun Lin, and Song Han.  Point-voxel CNN for efficient 3d deep learning.  _NeurIPS_ , 2019.
  * Loizou et al. [2020] Marios Loizou, Melinos Averkiou, and Evangelos Kalogerakis.  Learning part boundaries from 3d point clouds, 2020.
  * McInnes et al. [2018] Leland McInnes, John Healy, and James Melville.  Umap: Uniform manifold approximation and projection for dimension reduction, 2018.
  * Mildenhall et al. [2021] Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  _Commun. ACM_ , 65(1):99–106, 2021.
  * Mozaffar et al. [2021] Mojtaba Mozaffar, Shuheng Liao, Xiaoyu Xie, Sourav Saha, Chanwook Park, Jian Cao, Wing Kam Liu, and Zhengtao Gan.  Mechanistic artificial intelligence (mechanistic-AI) for modeling, design, and control of advanced manufacturing processes: Current state and perspectives.  _Journal of Materials Processing Technology_ , page 117485, 2021.
  * Powell [1964] Michael JD Powell.  An efficient method for finding the minimum of a function of several variables without calculating derivatives.  _The computer journal_ , 7(2):155–162, 1964.
  * Qi et al. [2017a] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  PointNet: Deep learning on point sets for 3d classification and segmentation.  In _CVPR_ , 2017a.
  * Qi et al. [2017b] Charles Ruizhongtai Qi, Li Yi, Hao Su, and Leonidas J Guibas.  PointNet++: Deep hierarchical feature learning on point sets in a metric space.  _NeurIPS_ , 2017b.
  * Rahaman et al. [2019] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville.  On the spectral bias of neural networks.  In _International Conference on Machine Learning_ , pages 5301–5310. PMLR, 2019.
  * Ranganathan [2004] Ananth Ranganathan.  The levenberg-marquardt algorithm.  _Tutoral on LM algorithm_ , 11(1):101–110, 2004\.
  * Rethage et al. [2018] Dario Rethage, Johanna Wald, Jurgen Sturm, Nassir Navab, and Federico Tombari.  Fully-convolutional point networks for large-scale point clouds.  In _ECCV_ , 2018.
  * Romanengo et al. [2023] Chiara Romanengo, Andrea Raffo, Silvia Biasotti, and Bianca Falcidieno.  Recognising geometric primitives in 3d point clouds of mechanical CAD objects.  _Computer-Aided Design_ , page 103479, 2023.
  * Roweis and Saul [2000] Sam T. Roweis and Lawrence K. Saul.  Nonlinear dimensionality reduction by locally linear embedding.  _Science_ , 290(5500):2323–2326, 2000.
  * Saporta and Sharf [2022] Tsahi Saporta and Andrei Sharf.  Unsupervised recursive deep fitting of 3d primitives to points.  _Computers & Graphics_, 102:289–299, 2022.
  * Sharma et al. [2020a] Gopal Sharma, Rishabh Goyal, Difan Liu, Evangelos Kalogerakis, and Subhransu Maji.  Neural shape parsers for constructive solid geometry.  _IEEE TPAMI_ , 44(5):2628–2640, 2020a.
  * Sharma et al. [2020b] Gopal Sharma, Difan Liu, Subhransu Maji, Evangelos Kalogerakis, Siddhartha Chaudhuri, and Radomír Měch.  ParSeNet: A parametric surface fitting network for 3d point clouds.  In _ECCV_ , 2020b.
  * Sitzmann et al. [2020] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  In _Advances in Neural Information Processing Systems_ , pages 7462–7473. Curran Associates, Inc., 2020.
  * Tenenbaum et al. [2000] Joshua B. Tenenbaum, Vin de Silva, and John C. Langford.  A global geometric framework for nonlinear dimensionality reduction.  _Science_ , 290(5500):2319–2323, 2000.
  * Uy et al. [2022] Mikaela Angelina Uy, Yen-Yu Chang, Minhyuk Sung, Purvi Goel, Joseph G Lambourne, Tolga Birdal, and Leonidas J Guibas.  Point2cyl: Reverse engineering 3d objects from point clouds to extrusion cylinders.  In _CVPR_ , 2022.
  * Vasu et al. [2021] Subeesh Vasu, Nicolas Talabot, Artem Lukoianov, Pierre Baqué, Jonathan Donier, and Pascal Fua.  HybridSDF: Combining deep implicit shapes and geometric primitives for 3d shape representation and manipulation.  _arXiv preprint arXiv:2109.10767_ , 2021.
  * Wang et al. [2012] Jun Wang, Dongxiao Gu, Zeyun Yu, Changbai Tan, and Laishui Zhou.  A framework for 3d model reconstruction in reverse engineering.  _Computers & Industrial Engineering_, 63(4):1189–1200, 2012.
  * Wang et al. [2019] Yue Wang, Yongbin Sun, Ziwei Liu, Sanjay E Sarma, Michael M Bronstein, and Justin M Solomon.  Dynamic graph cnn for learning on point clouds.  _ACM TOG_ , 38(5):1–12, 2019.
  * Willis et al. [2022] Karl DD Willis, Pradeep Kumar Jayaraman, Hang Chu, Yunsheng Tian, Yifei Li, Daniele Grandi, Aditya Sanghi, Linh Tran, Joseph G Lambourne, Armando Solar-Lezama, et al.  JoinABLe: Learning bottom-up assembly of parametric CAD joints.  In _CVPR_ , 2022.
  * Yan et al. [2021] Siming Yan, Zhenpei Yang, Chongyang Ma, Haibin Huang, Etienne Vouga, and Qixing Huang.  HPNet: Deep primitive segmentation using hybrid representations.  In _ICCV_ , 2021.
  * Yu et al. [2022] Fenggen Yu, Zhiqin Chen, Manyi Li, Aditya Sanghi, Hooman Shayani, Ali Mahdavi-Amiri, and Hao Zhang.  CAPRI-net: learning compact CAD shapes with adaptive primitive assembly.  In _CVPR_ , 2022.
  * Zhao et al. [2021] Hengshuang Zhao, Li Jiang, Jiaya Jia, Philip HS Torr, and Vladlen Koltun.  Point transformer.  In _ICCV_ , 2021.

Supplementary Material

##  6 Reproducibility

Please refer to the attached supplementary zip file for the source code of our
method. See README.md for instructions on setting up the environment and
running the code.

##  7 Notation

Throughout the paper, we use several terms and concepts as follows. A CAD
model represents a complex object in a format that permits editing in CAD
software. Boundary representation (B-REP) is one such representation; it
involves analytic surfaces, edges, and corners222 Related terms in literature:
“surfaces” := “faces”, “patches”; “edges” := “curves”, “contours”; “corners”
:= “points”, “vertices”, “endpoints”. , and explicit topological relationships
between them. Typically, such relationships take the form of “an edge is an
intersection of two surfaces” and “a corner is an intersection of two edges”,
stated as adjacency matrices. The manifolds of surfaces and edges have
intrinsic dimensionality of 2 and 1 respectively. Point cloud refers to the
unordered set of 3D points, obtained by scanning the surface of a real world
object or by sampling a CAD model for simulation purposes. Semantic classes
are the geometric types that we consider for representing CAD model surfaces.
These typically include a plane, sphere, cylinder, cone, and freeform (spline)
surfaces. Segmentation refers to assigning a semantic class to each point of
the point cloud. Clusters are sets of points belonging to the same semantic
class and instance; they group points belonging to distinct surfaces.
Clustering information can be a byproduct of segmentation, an algorithmic
stage, or given as ground truth. Primitives are inherently low-dimensional
parametric models representing analytical shapes of semantic classes, and that
can be used to fit clusters with least squares.

##  8 Primitive parameterization and fitting

Our method can handle several primitive types, including planes, spheres,
cylinders, and cones [24, 5].

#### Plane

In the context of 3D geometry, a plane can be represented by a vector
$\mathbf{n}\in\mathbb{R}^{3}$, which is the unit normal to the plane, and a
scalar $d$ that determines the distance of the plane from the origin.
Mathematically, a plane can be denoted by a tuple $(\mathbf{n},d)$, where
$\mathbf{n}$ satisfies $|\mathbf{n}|=1$. Any point
$\mathbf{p}\in\mathbb{R}^{3}$ lies on the plane if and only if
$\mathbf{n}^{\mathrm{T}}\mathbf{p}=d$.

The following steps can be employed to perform plane fitting to a collection
of 3D points:

  1. 1.

Calculate the centroid of the points:
$\mathbf{c}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{p_{i}}$;

  2. 2.

Compute the matrix:
$\mathbf{M}=\frac{1}{n}\sum_{i=1}^{n}(\mathbf{p_{i}}-\mathbf{c})(\mathbf{p_{i}}-\mathbf{c})^{T}$;

  3. 3.

Calculate the eigenvalues and eigenvectors of $\mathbf{M}$. Select the
eigenvector corresponding to the smallest eigenvalue as the normal vector
$\mathbf{n}$. Compute the distance d using the normal vector $\mathbf{n}$ and
the centroid c: $d=\mathbf{n}^{T}\mathbf{c}$.

#### Sphere

A sphere can be represented by a tuple $\left(\mathbf{c},r\right)$, where
$\mathbf{c}\in\mathbb{R}^{3}$ denotes the center of the sphere and
$r\in\mathbb{R}$ denotes its radius. Any point $\mathbf{p}\in\mathbb{R}^{3}$
lies on the sphere if and only if $||p-\mathbf{c}||=r$. The error function to
be minimized is as follows:

$\displaystyle
E\left(\mathbf{c},r\right)=\sum_{i=1}^{M}\left(\|\mathbf{p}_{i}-\mathbf{c}\|-r\right)^{2}.$
(3)

Solving $\frac{\partial E}{\partial r}=0$ for $r$ yields

$\displaystyle r=\frac{1}{M}\sum_{i=1}^{M}\|\mathbf{p}_{i}-\mathbf{c}\|.$ (4)

Solving $\frac{\partial E}{\partial\mathbf{c}}=0$ for $\mathbf{c}$ gives

$\displaystyle\mathbf{c}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{p}_{i}+r\cdot\frac{1}{M}\sum_{i=1}^{M}\frac{\partial(\|\mathbf{p}_{i}-\mathbf{c}\|)}{\partial\mathbf{c}}.$
(5)

A fixed point iteration can solve the equations and subsequently determine the
parameters $(\mathbf{c},r)$.

#### Cylinder

An infinite cylinder $(\mathbf{a},\mathbf{c},r)$ is characterized by a point
$\mathbf{c}$, a unit-length direction vector $\mathbf{a}$ that defines its
axis and a radius $r$. Any point $\mathbf{p}\in\mathbb{R}^{3}$ lies on the
cylinder if and only if
$\left(\mathbf{p}-\mathbf{c}\right)^{T}\left(\mathbf{I}-\mathbf{a}\mathbf{a}^{T}\right)\left(\mathbf{p}-\mathbf{c}\right)=r$
The error function to be minimized is

$\displaystyle
E\left(\mathbf{a},\mathbf{c},r\right)\\!=\\!\sum_{i=1}^{M}\\!\left(\left(\mathbf{p}_{i}\\!-\\!\mathbf{c}\right)^{T}\\!\\!\left(\mathbf{I}\\!-\\!\mathbf{a}\mathbf{a}^{T}\right)\\!\left(\mathbf{p}_{i}\\!-\\!\mathbf{c}\right)-r^{2}\right)^{2},$
(6)

where $\mathbf{I}$ is the identity matrix. Solving $\frac{\partial E}{\partial
r^{2}}=0$ leads to

$\displaystyle
r=\left(\frac{1}{M}\sum_{i=1}^{M}(\mathbf{c}-\mathbf{p}_{i})^{T}(\mathbf{I}-\mathbf{a}\mathbf{a}^{T})(\mathbf{c}-\mathbf{p}_{i})\right)^{\frac{1}{2}}.$
(7)

By solving $\frac{\partial E}{\partial\mathbf{c}}=0$, we can obtain an
equation for $\mathbf{c}$:

$\displaystyle\mathbf{c}=\frac{\tilde{\mathbf{A}}^{-1}\mathbf{A}}{\operatorname{tr}(\hat{\mathbf{A}}\mathbf{A})}\left(\frac{1}{M}\sum_{i=1}^{M}\left(\mathbf{p}_{i}^{T}\tilde{\mathbf{A}}\mathbf{p}_{i}\right)\mathbf{p}_{i}\right),$
(8)

where

$\displaystyle\mathbf{A}=\tilde{\mathbf{A}}\left(\frac{1}{M}\sum_{i=1}^{M}\mathbf{p}_{i}\mathbf{p}_{i}^{T}\right)\tilde{\mathbf{A}},$
(9)

and $\hat{\mathbf{A}}=\mathbf{S}\mathbf{A}\mathbf{S}^{T}$,
$\tilde{\mathbf{A}}=\mathbf{I}-\mathbf{a}\mathbf{a}^{T}$, $\mathbf{S}$ is the
skew symmetric matrix of $\mathbf{a}$. Putting them back into the error
function yields:

$\displaystyle G\left(\mathbf{a}\right)$
$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\biggl{[}\mathbf{p}_{i}^{T}\tilde{\mathbf{A}}\mathbf{p}_{i}-\frac{1}{M}\sum_{j=1}^{M}\mathbf{p}_{j}^{T}\tilde{\mathbf{A}}\mathbf{p}_{j}\,-$
(10)
$\displaystyle-\frac{2\mathbf{p}_{i}^{T}\hat{\mathbf{A}}}{\operatorname{Tr}\left(\hat{\mathbf{A}}\mathbf{A}\right)}\left(\frac{1}{M}\sum_{j=1}^{M}\left(\mathbf{p}_{j}^{T}\tilde{\mathbf{A}}\mathbf{p}_{j}\right)\mathbf{p}_{j}\right)\biggr{]}^{2}.$

A Powell optimizer [32] is then employed to locate the global minimum of the
function $G(\mathbf{a})$. Once $\mathbf{a}$ has been determined, the center
$\mathbf{c}$ and radius $r$ of the corresponding circle can be easily obtained
through Eq. 7 and 8.

#### Cone

We parameterize an infinite cone using the set of parameters
$(\mathbf{v},\mathbf{a},\theta)$, where $\mathbf{v}\in\mathbb{R}^{3}$ denotes
the apex point, $\mathbf{a}\in\mathbb{R}^{3}$ denotes a unit axis direction
vector, and $\theta\in(0,\pi/2)$ represents half the angle of the cone. Any
point $\mathbf{p}$ that lies on it satisfies
$\mathbf{a}\cdot\frac{\mathbf{p}-\mathbf{v}}{||\mathbf{p}-\mathbf{v}||}=\cos(\theta)$,
which can be written in a quadratic form as
$(\mathbf{p}-\mathbf{v})^{T}(\cos(\theta)^{2}\mathbf{I}-\mathbf{a}\mathbf{a}^{T})(\mathbf{p}-\mathbf{v})=0$.
Thus, the error function can be defined as follows:

$\displaystyle\\!\\!E\\!\left(\mathbf{v},\\!\mathbf{a},\\!\theta\right)\\!=\\!\sum_{i=1}^{M}\\!\left(\\!\left(\mathbf{p}_{i}\\!-\\!\mathbf{v}\right)^{T}\\!\\!\left(\\!\cos\left(\theta\right)^{2}\\!\mathbf{I}\\!-\\!\mathbf{a}\mathbf{a}^{T}\\!\right)\\!\left(\mathbf{p}_{i}\\!-\\!\mathbf{v}\right)\\!\right)^{2}\\!\\!.$
(11)

This least-square problem can be solved efficiently with the Levenberg-
Marquardt [36] method.

##  9 Algorithm of topological reconstruction

We present a pseudocode-based methodology for reconstructing topological
structures derived from $M$ surfaces fitted with potentially infinite
primitives, forming a set $\\{\mathbf{S}^{0}_{k}\\}_{M}$. We denote the stage
of surface processing with superscript. We trim each primitive to form a
margin of width $\epsilon$ around the input points, then employ tessellation
and triangulation meshing algorithm. As a result, we obtain a set of finite
extended surfaces denoted as $\\{\mathbf{S}^{1}_{k}\\}_{M}$. We generate poly-
line edges $\\{\mathbf{E}^{1}_{k}\\}_{K}$ by identifying intersecting surfaces
and computing pairwise intersections. And we trim the surfaces by the edges,
thus obtaining $\\{\mathbf{S}^{2}_{k}\\}_{M}$. Similarly, we intersect
adjacent poly-line edges to obtain corner points $\\{\mathbf{C}_{k}\\}_{L}$,
subsequently trimming edges accordingly to get final edges
$\\{\mathbf{E}^{2}_{k}\\}_{K}$. See Alg. 1.

Algorithm 1 Point2CAD model reconstruction

1:Input: $M$ surfaces $\\{\mathbf{S}^{0}_{k}\\}_{M}$ obtained from fitting
primitives to point clusters $\\{\mathbf{P}_{k}\\}_{M}$

2:Output: a CAD model with $M$ surfaces $\\{\mathbf{S}^{2}_{k}\\}_{M}$, $K$
edges $\\{\mathbf{E}^{2}_{k}\\}_{K}$ and $L$ corners
$\\{\mathbf{C}_{k}\\}_{L}$

3:for $i\in 1..M$ do

4: Trim $\mathbf{S}^{0}_{i}$ by $\epsilon$ to input points:
$\mathbf{S}^{1}_{i}=\operatorname{Trim}(\mathbf{S}^{0}_{i}|\mathbf{P}_{i},\epsilon)$

5:end for

6:for $i\in 1..M$ do

7: get edges on $\mathbf{S}^{1}_{i}$:
$\\{\mathbf{E}^{1}_{i,r}\\}_{R_{i}}=\\{\mathbf{S}^{1}_{j}\\}_{j\neq
i}\cap\mathbf{S}^{1}_{i}$

8: trim $\mathbf{S}^{1}_{i}$ by edges:
$\mathbf{S}^{2}_{i}=\operatorname{Trim}(\mathbf{S}^{1}_{i}|\\{\mathbf{E}^{1}_{i,r}\\}_{R_{i}})$

9:end for

10:for each pair $(\mathbf{E}_{p},\mathbf{E}_{q})$ in intersection edges do

11: obtain the corners $\mathbf{C}_{pq}=\mathbf{E}_{p}\cap\mathbf{E}_{q}$

12:end for

13:for $i\in 1..L$ do

14: if any $\mathbf{C}_{s}\in\mathbf{E}_{i}$ then

15: trim it by the corners:
$\mathbf{E}^{2}_{i}=\operatorname{Trim}(\mathbf{E}^{1}_{i}|\mathbf{C}_{s})$

16: end if

17:end for

Two distinct cases exist of using the $\operatorname{Trim}$ operation on
surfaces. In the first case, trimming is performed based on a distance
threshold, retaining only the portion of the infinite primitives near the
input points. In the second case, we trim the surfaces by considering the
intersection of their triangle-mesh representations with edges. We employ
connected component analysis, whereby a pair of faces is considered connected
if a path exists between them that does not cross an edge obtained through the
surface intersection. We then discard whole connected components based on the
distance of their members to the original point cloud. The
$\operatorname{Trim}$ operation on edges involves a similar subdivision of the
edge into connected segments by corners and retaining segments close to the
input points.


# Edge-sharing quasi-one-dimensional cuprate fragments in optimally
substituted Cu/Pb apatite

Katherine Inzani1 1School of Chemistry, University of Nottingham, University
Park, Nottingham NG7 2RD, United Kingdom  John Vinson2 2Material Measurement
Laboratory, National Institute of Standards and Technology, Gaithersburg,
Maryland, 20899, USA  Sinéad M. Griffin3,4 3Materials Sciences Division,
Lawrence Berkeley National Laboratory, Berkeley, California, 94720, USA
4Molecular Foundry Division, Lawrence Berkeley National Laboratory, Berkeley,
California, 94720, USA

(August 28, 2024)

###### Abstract

The flurry of theoretical and experimental studies following the report of
room-temperature superconductivity at ambient pressure in Cu-substituted lead
apatite CuxPb10-x(PO4)6O (‘LK99’) have explored whether and how this system
might host strongly correlated physics including superconductivity. While
first-principles calculations at low doping ($x\approx 1$) have indicated a
Cu-$d^{9}$ configuration coordinated with oxygen giving rise to isolated,
correlated bands, its other structural, electronic, and magnetic properties
diverge significantly from those of other known cuprate systems. Here we find
that higher densities of ordered Cu substitutions can result in the formation
of contiguous edge-sharing Cu-O chains, akin to those found in some members of
the cuprate superconductor family. Interestingly, while such quasi-one-
dimensional edge-sharing chains are typically ferromagnetically coupled along
the chain, we find an antiferromagnetic ground-state magnetic order for our
cuprate fragments which is in proximity to a ferromagnetic quantum critical
point. This is a result of the elongated Cu-Cu distance in Cu-substituted
apatite that leads to larger Cu-O-Cu angles supporting antiferromagnetism,
which we demonstrate to be controllable by strain. Finally, our electronic
structure calculations confirm the low-dimensional nature of the system and
show that the bandwidth is driven by the Cu-O plaquette connectivity,
resulting in an intermediate correlated regime.

††preprint: APS/123-QED

##  I Introduction

Cuprate systems exhibit a rich playground of strongly-correlated physics
encompassing phenomena from high-temperature superconductivity to frustrated,
low dimensional quantum spin systems [1, 2]. Despite this diversity of
emergent phenomena, several key features are shared across these cuprate
systems. The stoichiometric compounds are characterized as Mott or charge-
transfer insulators with $S=\frac{1}{2}$ moments on Cu ions stemming from a
Cu-$d^{9}$ configuration. These Cu ions are coordinated with oxygen to form
CuO4 square planes that are typically antiferromagnetically coupled. Electron-
or hole-doping these antiferromagnetic insulators can suppress the magnetic
order and result in superconductivity. It is these superexchange-driven
correlations that are suggested to give rise to potential $d$-wave
superconductivity [3, 4].

The fundamental building blocks of the cuprates are square planar coordinated
CuO4 that form either corner- or edge-sharing networks. Corner-sharing chains,
found in cuprate superconductors and spin ladders, exhibit strong nearest
neighbour (NN) exchange interactions (100 meV to 160 meV) facilitated through
an $\approx 180^{\circ}$ Cu-O-Cu pathway. Conversely, edge-sharing chains have
weaker NN exchange interactions due to their $\approx 90^{\circ}$ Cu-O-Cu
pathway, resulting in a range of reported ferromagnetic and antiferromagnetic
behaviors dependent on the specific Cu-O-Cu bond angle. Moreover, reducing the
dimensionality of these $S=\frac{1}{2}$ systems, from the well-known 2D
antiferromagnetic Heisenberg model to 1D, results in a significant enhancement
of quantum fluctuations. In fact, Bethe predicted such strong quantum
fluctuations would prevent long-range magnetic order from forming in the 1D
$S=\frac{1}{2}$ antiferromagnetic Heisenberg chain [5]. The diversity in
structural configurations of cuprates, influenced by factors such as
connectivity, orbital overlap, and dimensionality, highlights the complexity
of the physical landscape inherent in cuprates, leading to a wide spectrum of
observed properties in these compounds.

Figure 1: Crystal structures with the primitive unit cells of the non-magnetic
cell shown by the box. The arrows depict the calculated ground-state magnetic
order. (a) Cu2Pb${}^{\mathrm{I}}_{2}$Pb${}^{\mathrm{II}}_{6}$(PO4)6O with
(half) two PbI sites substituted with Cu. (b)
Cu4Pb${}^{\mathrm{I}}_{0}$Pb${}^{\mathrm{II}}_{6}$(PO4)6O with all four PbI
sites substituted with Cu. (c) Orthorhombic Li2CuO2.

The recent report of possible room temperature superconductivity at ambient
pressures in Cu-substituted Pb apatite (also known as ‘LK99’) [6, 7] has
focused attention on their reported CuPb9(PO4)6O structure. Initial
calculations using Density Functional Theory (DFT) showed that Cu substitution
on the trigonal-prismatically coordinated Pb site results in a $d^{9}$
electronic configuration with a narrow, isolated, half-filled Cu-$d$ band in
its paramagnetic phase [8]. That work, and subsequent theoretical studies
predicted it to have a frustrated antiferromagnetic ground state at the DFT
level [9], however with very weak coupling. Appropriately incorporating strong
correlations through Dynamical Mean-Field Theory, showed it to lie in an
ultra-correlated regime ($U/W\approx 30-50$) which discounted known mechanisms
for electronically-driven superconductivity [10]. This motivates us to
consider how the bandwidth can be increased to a more moderately correlated
regime via chemical and physical modifications to the Cu-substituted apatite
framework. In this paper, we perform DFT calculations of Cu-substituted lead
apatite in higher Cu concentrations, namely CuxPb10-x(PO4)6O (for $x=2,4$),
identifying edge-shared square planar arrays of CuO4 plaquettes. We examine
the magnetic and electronic phase diagram of this system, comparing its key
features with the prototypical examples of cuprates with (i) edge-sharing
motifs, Li2CuO2 and (ii) corner-sharing motifs, La2CuO4.

##  II Methods

All of our DFT calculations used the Vienna Ab initio Simulation Package
(VASP) [11, 12, 13, 14], with further details given in the SI. To account for
Cu-$d$ localization, we applied a Hubbard-U correction of 8 eV which was
benchmarked in recent work with meta-GGA functionals and calculation of self-
consistent U [15]. Electronic structure plots were generated using the sumo
software package [16], and Fermi surfaces were generated using iFermi [17]. We
also benchmarked our calculations by repeating the electronic structure
calculation for a representative system with the meta-GGA r2SCAN functional
[18], finding qualitatively similar results (see SI).

##  III Results

###  III.1 Structural Properties

We consider here the family of apatites given by
CuxPb${}^{\mathrm{I}}_{4-x}$Pb${}^{\mathrm{II}}_{6}$(PO4-y)6X2, where we limit
the Cu substitution to the PbI site and X = O2-. Proper charge is maintained
by $y=1/6$, reflecting evidence that two O2- ions in the channel are balanced
with O vacancies elsewhere [19]. Our choice of only substituting Cu on site I
is motivated by the markedly different behavior of the Cu-$d$ states on site I
compared to site II. At $x=1$ substitution and ambient pressure synthesis,
neither site I and II occupancies are strongly favored [19]. However, Ca-doped
Pb apatite shows a strong affinity for ordered substitution on site I when
Ca:Pb ratio is 2:3 due to the mismatch in cation size, known as the hedyphane
sub-group of apatites, facilitating a smaller-than-trend unit cell compared to
the otherwise solid-solution behavior [20, 21, 22]. Theoretical investigations
using DFT support this approach, indicating that site I is preferred for $x=2$
substitution (Cu2Pb8(PO4)6) [23], and, for $x=1$ substitution, that site I is
preferred at high pressures [15, 24]. We therefore calculate Cu substitution
on site I for both $x=2$ ‘Cu-2’ and $x=4$ ‘Cu-4’, as depicted in Fig. 1,
noting that the original papers considered $x=1$.

The optimized structures for Cu-2 and Cu-4 are shown in Fig. 1 where both
include an oxygen vacancy (see SI for structure files). We note that including
oxygen vacancies seems crucial for the formation of Cu-O chains – for
calculations without oxygen vacancies and a single channel oxygen our
calculations find that contiguous plaquettes do not form (see SI). This
suggests that the formation of continguous plaquettes is sensitive to Cu site-
selectivity, Cu clustering, and oxygen vacancies. We find a significant
lattice collapse for both cases as a result of the substitution of the smaller
Cu2+ ion (0.73 Å) for Pb2+ site (1.19 Å) – the out-of-plane lattice parameters
reduce to $c=6.95$ Å and $c=5.99$ Å for Cu-2 and Cu-4, respectively. Cu
substitution on the PbI site distorts the original trigonal antiprism oxygen
coordination to now being square planar, resulting in a new monoclinic space
group of Pm (#6). Looking closely at the new Cu coordination, we find that the
$c$-lattice reduction facilitates the formation of contiguous edge-sharing
CuO4 plaquettes. Such edge-sharing cuprate fragments are less commonly found
than their corner-sharing counterparts, and are typically seen in high-
pressure phases such as SrCuO2 and in orthorhombic Li2CuO2 (Fig. 1(c)).

The reduction in $c$ lattice length is likely necessary to support the
formation of edge-sharing CuO4 plaquettes. Known edge-sharing cuprates have a
Cu-Cu separation that ranges from 2.86 Å in Li2CuO2 to 2.94 Å in CuGeO3, as
compared to our results of 3.48 Å in Cu-2 and 2.97 Å in Cu-4. We find that the
Cu-O bond distances range from 2.07 Å to 2.21 Å in Cu-2 to 1.95 Å to 2.01 Å in
Cu-4, compared to 1.96 Å in LiCuO2 or 1.93 Å in CuGeO3.

This variation in Cu-O bond lengths of the apatites is due in part to the fact
that the CuO4 plaquettes in Cu-2 and Cu-4 are puckered out-of-phase where the
coordinating oxygen are staggered in the plane perpendicular to the chain (see
Fig. 4). This puckering is also observed in some high-pressure phases of
cuprates such as incommensurate Sr0.73CuO2 [25].

Figure 2: Calculated non-spin-polarized electronic band structures for (a)
Li2CuO2, and (b) La2CuO4, (c) Cu-2 Cu2Pb8(PO4)6O, and (d) Cu-4 Cu4Pb6(PO4)6O.
The red and blue show the orbital projections of the Cu-$d$ and O-$p$ states
respectively. The Fermi level is set to 0 eV in all plots and is marked by the
dashed line. Note that the energy range for (b) is greater than the other
three – for comparison the energy range of (a), (c) and (d) is marked by the
dotted lines. $E_{\textrm{F}}\pm 0.1$ eV is also marked by dotted lines on (c)
and (d).

###  III.2 Electronic Structure

We first present the calculated electronic structures without including spin-
polarization which is more relevant for temperatures above the magnetic
Curie/Néel temperature (spin-polarized calculations are given in the SI).
These are plotted in Fig. 2 for the two Cu-ordered cases Cu-2 and Cu-4, and
for the known cuprates Li2CuO2 and La2CuO4 for a region around the Fermi
level. As expected from Cu-$d^{9}$ in a square-planar crystal field
environment in all structures, we find half-filled bands at the Fermi level
dominated by the highest-energy Cu-$d$ orbital. We note that the nonmagnetic
primitive unit cell of Li2CuO2 contains only one Cu ion whereas that of Cu-2
and Cu-4 contain two and four respectively, resulting in the half-filling of
two and four bands at the Fermi level.

However, beyond this, the electronic structures of these cuprate compounds
diverge. For Li2CuO2, the splitting of the hole orbital is much greater than
that of La2CuO4, Cu-2, or Cu-4 resulting in an isolated band separated by over
1 eV from the rest of the valence manifold. This is in contrast to the other
three compounds whose hole band hybridizes with the top of the valence band.
The bandwidth of Li2CuO2 is $\approx$1.1 eV, as a result of hopping along
edge-sharing Cu-O chains with a Cu-Cu separation of 2.83 Å. In contrast, the
Cu-Cu separation for Cu-2 and Cu-4 are 3.47 Å and 2.97 Å resulting in reduced
bandwidths of approximately 0.5 eV and 1 eV, respectively. These bandwidths
are greatly increased from the $x=1$ case reported previously of around 0.1 eV
[8], pushing these heavily substituted Cu-2 and Cu-4 apatites into a more
intermediately correlated regime where $U/W\approx 1-10$ as desired. The
corner-shared La2CuO4 has a bandwidth of over 3 eV, which at first glance
might be surprising given its Cu-Cu distance of 3.95 Å, the largest of all of
the cases considered here. The difference is in the Cu-O-Cu bond angle. For
the case of the previously discussed edge-shared Cu-O plaquettes, the Cu-$d$
and O-$p$ are orthogonal with a Cu-O-Cu bond angle of 90∘, resulting in a
strongly suppressed hopping. However, the corner-shared systems like La2CuO4
have significant Cu-$d_{x^{2}-y^{2}}$ hybridization with O-$p$ resulting in
increased bandwidth.

For corner-sharing square plaquettes, like those in La2CuO4, the orbital order
is dominated by crystal-field splitting of the square planar array which
results in the hole orbital being primarily Cu-$d_{x^{2}-y^{2}}$ with some
O-$p$ hybridization (see SI. Fig. S3(b)) as reported previously. However, for
the edge-sharing CuO2 plaquettes, metal-metal bonding results in a
redistribution of the orbitals such that the in-plane $\sigma$-type Cu-$d$
orbital is now half-filled. This corresponds to the orbital that is parallel
to the Cu-O-Cu chain, and perpendicular to the interchain direction [26]. For
the case of Li2CuO2, the hole band therefore has Cu-$d_{xz}$ character with
significant hybridization found with O-pz, as previously reported [27].
Similarly for both Cu-2 and Cu-4, the chain runs in the out-of-plane
direction, resulting in a hole band of Cu-$d_{xz}$ with some O-$p_{z}$ (see
SI. Fig. S3(c,d)).

In Fig. 3 we plot the non-spin-polarized Fermi surfaces for Cu-2 and Cu-4 for
a range of values of the Fermi level. Interestingly, for both cases we find
nested Fermi surfaces and a strong dependence of the dimension and nesting on
the position of the Fermi level. For instance, Cu-4 transitions from being 3D
to more 2D in going from 0.1 eV below the Fermi level to the Fermi level. The
Fermi surface character also changes from comprised concentric cylinders
surrounding the center of the BZ edge ($X$ point) at the Fermi level to two
separate cylinders centered on the BZ edge and corner at 0.1 eV above the
Fermi level. This sensitivity of the details of the Fermi surface with doping
also suggests strong potential spin frustration and fluctuations in these
systems.

Figure 3: Calculated non-spin-polarized Fermi surfaces for Cu-2 at (a)
$E_{\textrm{F}}-0.1$ eV, (b) $E_{\textrm{F}}$, and (c) $E_{\textrm{F}}+0.1$
eV, and for Cu-4 in (d), (e) and (f) respectively.

###  III.3 Magnetic Properties

Figure 4: Calculated energy difference between ferromagnetic and
antiferromagnetic coupling along the Cu-O chains for
Cu2Pb${}^{\mathrm{I}}_{2}$Pb${}^{\mathrm{II}}_{6}$(PO4)6O as a function of
average Cu-O-Cu bond angle when uniaxial out-of-plane strain is applied.
Negative values correspond to an antiferromagnetic ground state. The top panel
shows the evolution of apical oxygen bond lengths for a variety of
configurations corresponding to 94.6∘, 101.8∘ and 106.4∘.

In edge-sharing cuprates, the nature of frustrated interactions is deeply
intertwined with the linkage type of CuO4 plaquettes. With a Cu-O-Cu angle
adjusting from 180∘ in corner sharing to 90∘ in edge sharing, the nature of
the exchange interaction dramatically changes. According to the Goodenough-
Kanamori-Anderson rule, the Cu-Cu nearest-neighbor (NN) exchange interaction
is strong and antiferromagnetic across a 180∘ Cu-O-Cu bond, while a 90∘ bond
results in a weaker and ferromagnetic exchange. This implies that in edge-
sharing complexes, the antiferromagnetic superexchange interaction is
diminished due to the orthogonality of Cu-$d$ and O-$p$ orbitals. However, an
effective exchange interaction is still plausible through the Cu-O-O-Cu path
between the next-nearest neighbor (NNN) Cu ions, with the contribution of the
NNN exchange interaction being antiferromagnetic. This significant influence
of the NNN leads to magnetic frustration in edge-sharing systems. This then
results in diverse ground states dependent on the ratios of these competing
exchange interactions and the corresponding connection to the CuO4
connectivity.

We calculated the ground-state magnetic order of the Cu-4 structure in the
collinear limit finding it to be G-type with antiferromagnetic interactions
both along the chains and between the chains. We find this to be 0.9 meV/Cu
(5.3 meV/Cu) more favorable than ferromagnetic interactions along the chain
(and antiferromagnetic between chains), and 1.1 meV/Cu (5.0 meV/Cu) more
favorable than antiferromagnetic along chains that are ferromagnetically
coupled, indicating strong frustration. The values in parenthesis are those
calculated using the r2SCAN functional, which shows an enhancement of exchange
interactions, but maintains the ground-state AFM order. This AFM coupling
along the chain can be rationalized by the Goodenough-Kannamori rules that
states that in going from a Cu-O-Cu bond angle of 180∘ to 90∘, the exchange
interaction changes from antiferromagnetic to ferromagnetic. In our case the
Cu-O-Cu bond angles vary between 96∘ and 102∘, consistent with our DFT-
calculated antiferromagnetic spin chain order.

For Cu-2, we also find a G-type AFM ground state with AFM coupling both along
the chain and between chains. This is only 8 $\mu$eV/Cu (21 $\mu$eV/Cu for
r2SCAN) more favorable than AFM coupling along the chain and FM coupling
between the chains, as expected from the larger chain separation. For the case
of FM coupling along chains that are AFM coupled to each other, we find it to
be 7.1 meV/Cu (60 meV/Cu for r2SCAN) higher in energy than both AFM coupled.
Like the Cu-4 case, we find AFM exchange along the chains, independent of
functional, however for Cu-2, AFM further stabilized. This is due to the
larger $c$ lattice for Cu-2 which results in much larger Cu-O-Cu bond angles
ranging from 106∘ and 110∘. Given the strong dependence of exchange
interactions on the out-of-plane lattice parameter, $c$, which in turn
modifies the Cu-O-Cu bond angle, we next explicitly calculate how uniaxial
strain along $c$ can be used to control exchange along the chains. Given the
very weak interchain coupling for Cu-2, we neglect it for the strain
calculations, keeping it fixed as FM. In Fig. 4 we plot the calculated energy
difference between AFM and FM intrachain coupling for Cu-2 as a function of
the smallest Cu-O-Cu bond angle resulting from adjusting the $c$-lattice
parameter (parallel to the chain). The fully relaxed Cu-2 structure
corresponds to the right-most datapoint of 106.4∘. We confirm the strong
dependence of exchange interactions on bond angle, uncovering a transition
from AFM to FM via a quantum critical point.

We next compare these results with other reported edge-sharing cuprate
systems, in particular for Li2CuO2 where we have performed additional DFT
calculations for direct comparison. We confirm that Li2CuO2 has an overall
antiferromagnetic ground state, comprising ferromagnetic interactions along
the Cu-O chains with the chains then antiferromagnetically coupled to each
other. This is consistent with previous DFT calculations and neutron
diffraction experiments [28, 27] where we also find a significant fraction
($\approx$20%) of the ordered moment to be localized on the O-$p$ states. From
our structural optimization we find that the Cu-O-Cu bond angle is 94∘,
putting it closer to the ferromagnetic exchange predictions from the
Goodenough-Kannomori rules. Other systems with edge-sharing spin chains also
follow this trend where the Cu-O-Cu bond angle dictates the dominant exchange
interactions – this is typified in the example of high-pressure Sr0.73CuO2. It
adopts an incommensurate structure owing to the mismatch between the Sr ions
and Cu-O plaquette periodicities – this results in a distorted wave of the
Cu-O plaquettes with distances ranging from 1.85 Å to 2.05 Å. Like Li2CuO2, it
is also ferromagnetic along the Cu-O chains, which are antiferromagnetically
coupled to each other – the coupling strength along the chain can be modified
with the Ca substitution on the Sr site [29]. On the other hand, the corner-
sharing undoped La2CuO4 adopts a G-type antiferromagnetic order which we
calculate to be 68 meV/Cu and 158 meV/Cu more favorable than striped AFM and
FM orders respectively, consistent with its 180∘ Cu-O-Cu bond angle.

##  IV Discussion

The Cu-substituted apatites provide a new system for understanding one-
dimensional cuprate chains in a previously inaccessible coupling regime. The
large lattice vector (Cu-Cu separation) compared with other edge-sharing
cuprates leads to a larger Cu-O-Cu angle. In turn, this enabled a stronger
effective Cu-Cu nearest-neighbor interaction, stabilizing antiferromagnetic
ordering, in contrast to typically ferromagnetic edge-sharing cuprates. Our
calculations show that the Cu-4 structure is very close to a ferromagnetic
ground state, indicative of strongly frustrated interactions and proximity to
a quantum critical point. As shown in Fig. 4, this transition is dependent on
the Cu-O-Cu angle, and can be tuned by structural modifications.

The increased band width of the Cu-$d_{xy}$ states at the Fermi level seen in
Cu-2 and Cu-4 over the previously examined $x=1$ structure is an expected
result of halving the Cu-Cu separation along the $c$ direction. While the Cu-
Cu separation is large compared to other edge-sharing cuprates, which could be
expected to reduce the band width, we find the structures to be in the
intermediate coupling regime, favorable for superconductivity. This is in part
due to the distortion to larger Cu-O-Cu bond angles which is also seen to
stabilize an antiferromagnetic ground state. The oxygen atoms that constitute
the plaquettes are also bound to the PO4 units, the size of which, along with
the $c$ lattice constrains the Cu-O bond lengths and the Cu-O-Cu angle.
Substitutions on the B sites may allow for a further tuning of the correlation
and antiferromagnetic ordering.

The proposed structure requires substantially higher doping than has been
reported experimentally [6, 7, 19] and also with a high degree of site
selectivity to facilitate macroscopic cuprate chains. This is not
insurmountable as both solid-solution and ordered doping of apatite structures
are common. Reported synthesis attempts have been carried out under ambient
pressure and using oxygen as the channel species X. Recent theoretical work
has identified that pressure, strain, and channel species call all have a
strong effect on the site selectivity of the copper dopants [15, 24, 23].

#### Acknowledgements —

We are grateful to Alexander Balatsky and Masatoshi Imada for useful feedback.
This work was funded in part by the U.S. Department of Energy, Office of
Science, Office of Basic Energy Sciences, Materials Sciences and Engineering
Division under Contract No. DE-AC02-05-CH11231 within the Theory of Materials
program. Computational resources were provided by the National Energy Research
Scientific Computing Center and the Molecular Foundry, DOE Office of Science
User Facilities supported by the Office of Science, U.S. Department of Energy
under Contract No. DEAC02-05CH11231. The work performed at the Molecular
Foundry was supported by the Office of Science, Office of Basic Energy
Sciences, of the U.S. Department of Energy under the same contract. K.I.
acknowledges support from the EPSRC (EP/W028131/1).

Certain software is identified in in this paper in order to specify the
procedure adequately. Such identification is not intended to imply
recommendation or endorsement of any product or service by NIST, nor is it
intended to imply that it is necessarily the best available for the purpose.

## References

  * Bednorz and Müller [1986] J. G. Bednorz and K. A. Müller, Zeitschrift für Physik B Condensed Matter 64, 189 (1986).
  * Keimer et al. [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Scalapino et al. [1986] D. Scalapino, E. Loh Jr, and J. Hirsch, Phys. Rev. B 34, 8190 (1986).
  * Anderson [1987] P. W. Anderson, Science 235, 1196 (1987).
  * Bethe [1931] H. Bethe, Zeitschrift für Physik 71, 205 (1931).
  * Lee et al. [2023a] S. Lee, J.-H. Kim, and Y.-W. Kwon, arXiv preprint arXiv:2307.12008 (2023a).
  * Lee et al. [2023b] S. Lee, J. Kim, H.-T. Kim, S. Im, S. An, and K. H. Auh, arXiv preprint arXiv:2307.12037 (2023b).
  * Griffin [2023] S. M. Griffin, arXiv preprint arXiv:2307.16892 (2023).
  * Kurleto et al. [2023] R. Kurleto, S. Lany, D. Pashov, S. Acharya, M. van Schilfgaarde, and D. S. Dessau, arXiv preprint arXiv:2308.00698 (2023).
  * Yue et al. [2023] C. Yue, V. Christiansson, and P. Werner, arXiv preprint arXiv:2308.04976 (2023).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Comp. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Inzani et al. [2024] K. Inzani, J. Vinson, T. Bayaraa, J. Broad, N. Dale, B. Field, T. Huynh, G. Moore, and S. M. Griffin, In Preparation (2024).
  * Ganose et al. [2018] A. M. Ganose, A. J. Jackson, and D. O. Scanlon, Journal of Open Source Software 3, 717 (2018).
  * Ganose et al. [2021] A. M. Ganose, A. Searle, A. Jain, and S. M. Griffin, Journal of Open Source Software 6, 3089 (2021).
  * Furness et al. [2020] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, The Journal of Physical Chemistry Letters 11, 8208 (2020).
  * Puphal et al. [2023] P. Puphal, M. Y. P. Akbar, M. Hepting, E. Goering, M. I. A. A, and B. K. Nugroho, arXiv preprint arXiv:2308.06256 (2023).
  * Rouse et al. [1984] R. C. Rouse, P. J. Dunn, and D. R. Peacor, American Mineralogist 69, 920 (1984).
  * Bigi et al. [1989] A. Bigi, A. Ripamonti, S. Brückner, M. Gazzano, N. Roveri, and S. A. Thomas, Acta Crystallographica Section B 45, 247 (1989).
  * Bigi et al. [1991] A. Bigi, M. Gandolfi, M. Gazzano, A. Ripamonti, N. Roveri, and S. A. Thomas, J. Chem. Soc., Dalton Trans. pp. 2883–2886 (1991).
  * Shen et al. [2023] J. Shen, D. G. II, S. Shahabfar, Z. Li, D. Kang, S. Griesemer, A. Salgado-Casanova, T. chen Liu, C.-T. Chou, Y. Xia, et al., arXiv preprint arXiv:2308.07941 (2023).
  * Ogawa et al. [2023] K. Ogawa, K. Tolborg, and A. Walsh, ACS Energy Letters 8, 3941 (2023).
  * Schwer et al. [1997] H. Schwer, K. Conder, E. Kopnin, G. Meijer, R. Molinski, and J. Karpinski, Physica C: Superconductivity 282, 901 (1997).
  * Huang et al. [2011] H.-Y. Huang, N. A. Bogdanov, L. Siurakshina, P. Fulde, J. van den Brink, and L. Hozoi, Phys. Rev. B 84, 235125 (2011).
  * Weht and Pickett [1998] R. Weht and W. E. Pickett, Phys. Rev. Lett. 81, 2502 (1998).
  * Sapina et al. [1990] F. Sapina, J. Rodriguez-Carvajal, M. Sanchis, R. Ibanez, A. Beltran, and D. Beltran, Solid State Communications 74, 779 (1990).
  * Wu et al. [1999] H. Wu, Q.-Q. Zheng, X.-G. Gong, and H. Lin, Journal of Physics: Condensed Matter 11, 4637 (1999).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew et al. [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
  * Dudarev et al. [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).

##  V Supplemental Material

###  V.1 Calculation Details

We used the Vienna Ab initio Simulation Package (VASP) [11, 12, 13, 14] for
all of our Density Functional Theory (DFT) calculations with projector
augmented wave (PAW) pseudopotentials [30, 31]. We included the Pb 5d106s26p2,
Cu 3d104s1, P 3s23p5, O 2s22p4, Li 2s1, and La 5s25p65d16s2 as valence
electrons. For all calculations we used a plane-wave cut-off energy of 700 eV.
The Brillouin Zone (BZ) was sampled with a $6\times 6\times 8$ Gamma-centered
k-point grid for structural optimizations and a $10\times 10\times 12$ Gamma-
centered k-point grid for the density of states for the Cu-2 and Cu-4 systems.
For Li2CuO2 and La2CuO4, we used a $12\times 12\times 4$ Monkhorst-Pack
k-point grid. We used the generalized gradient approximation (GGA) based
exchange-correlation functional PBEsol [32], and applied a Hubbard-U of 8 eV
to the Cu-$d$ states within the Dudarev approach [33], unless otherwise
stated. Meta-GGA functional calculations used the r2SCAN functional [18]. The
electronic convergence criterion is set to $10^{-6}$ eV and the force
convergence criterion is set to 0.01 eV / Å.

###  V.2 Choice of Exchange-Correlation Functional

We benchmarked our choice of exchange-correlation functional by comparing the
DFT+U (U = 8 eV) magnetic and electronic structure properties with a select
number of meta-GGA calculations using the r2SCAN functional. We find
qualitatively similar results for the magnetic ground states regardless of
functional choice, however the calculated exchange coupling strengths are
greater with r2SCAN (see discussion in main text). We also calculate the
electronic band structure for Cu-2 with r2SCAN finding, again, good
qualitative agreement with very similar hole and electron pockets near the
Fermi level. However, we find small shifts in the relative band filling at the
Fermi level – for instance the O-$p$ states are above $E_{F}$ with DFT+U, see
Fig. 2(c), while these are below $E_{F}$ with r2SCAN (Fig. S1). This suggests
that correctly treating the O-$p$ states will also be important for accurately
describing the low-energy physics of these systems.

Figure S1: Non-spin-polarized electronic band structure of Cu-2 calculated
with r2SCAN on PBEsol geometry. The Fermi level is set to 0 eV and is marked
by the dashed line. Figure S2: Calculated spin-polarized electronic band
structure (left) and corresponding density of states (right) for (a)
Orthorhombic Li2CuO2, and (b) Cu4Pb6(PO4)6(O), ‘Cu-4’ where the Cu are all on
the Pb-1 sites as depicted in Fig.1. The spin-up bands are depicted in solid
orange, and the spin-down bands are dashed blue. The total density of states
is shaded grey with projections shown of the Cu-$d$ orbitals (pink) and its
neighboring O-$p$ orbitals (green). In both plots the Fermi level is set to 0
eV and is marked by the dashed line.  Figure S3: Calculated non-spin-polarized
electronic band structures (left) and orbital-projected densities of states
(right) for (a) Li2CuO2, and (b) La2CuO4, (c) Cu-2 Cu2Pb8(PO4)6O, and (d) Cu-4
Cu4Pb6(PO4)6O. The Fermi level is set to 0 eV in all plots and is marked by
the dashed line.

##  VI Structure Files

The resulting relaxed structures in POSCAR format are given below. The
specific geometry is given in the headers.

    Cu-2 with oxygen vacancy; AFM
        1.000000000000000
         9.90522    0.73241    0.00000
        -6.06072    8.01068    0.00000
         0.00000    0.00000    6.95031
       Cu   Pb   P    O
         2     8     6    25
    Direct
      0.30643  0.71928  0.99956
      0.30643  0.71928  0.50044
      0.67644  0.32829  0.99414
      0.67644  0.32829  0.50586
      0.74168  0.70285  0.25000
      0.24732  0.26673  0.75000
      0.20552  0.99261  0.25000
      0.78750  0.99542  0.75000
      0.96639  0.19554  0.25000
      0.98062  0.77444  0.75000
      0.39661  0.40388  0.25000
      0.60204  0.62780  0.75000
      0.59361  0.01747  0.25000
      0.36831  0.95856  0.75000
      0.98601  0.57945  0.25000
      0.04285  0.41098  0.75000
      0.02251  0.01514  0.75000
      0.94150  0.93615  0.25000
      0.34740  0.51899  0.25000
      0.64208  0.49746  0.75000
      0.46666  0.83515  0.25000
      0.50027  0.13556  0.75000
      0.17308  0.65866  0.25000
      0.86526  0.36459  0.75000
      0.42389  0.56266  0.75000
      0.50784  0.10887  0.25000
      0.43392  0.84590  0.75000
      0.89636  0.39643  0.25000
      0.15998  0.59702  0.75000
      0.39417  0.31753  0.07069
      0.69093  0.74269  0.92464
      0.70577  0.07218  0.07152
      0.26185  0.90985  0.93319
      0.93494  0.63223  0.07137
      0.07539  0.34189  0.92949
      0.69093  0.74269  0.57536
      0.39417  0.31753  0.42931
      0.26185  0.90985  0.56681
      0.70577  0.07218  0.42848
      0.07539  0.34189  0.57051
      0.93494  0.63223  0.42863

    Cu-4 with oxygen vacancy; AFM
        1.000000000000000
        10.26014   -0.19929    0.00000
        -7.09771    9.12119    0.00000
         0.00000    0.00000    5.99455
       Cu   Pb   P    O
         4     6     6    25
      0.28115  0.71803  0.99746
      0.67207  0.25768  0.99870
      0.67207  0.25768  0.50130
      0.28115  0.71803  0.50254
      0.72766  0.67239  0.25000
      0.31283  0.32332  0.75000
      0.22296  0.03758  0.25000
      0.84480  0.01509  0.75000
      0.99048  0.21948  0.25000
      0.97660  0.79932  0.75000
      0.48079  0.37543  0.25000
      0.58274  0.67734  0.75000
      0.59368  0.01267  0.25000
      0.36291  0.96230  0.75000
      0.93906  0.55758  0.25000
      0.01974  0.42010  0.75000
      0.10394  0.06577  0.75000
      0.95491  0.87840  0.25000
      0.32658  0.36613  0.25000
      0.61026  0.56134  0.75000
      0.43884  0.84251  0.25000
      0.51620  0.13026  0.75000
      0.12783  0.61861  0.25000
      0.82971  0.35047  0.75000
      0.39966  0.61062  0.75000
      0.53425  0.11132  0.25000
      0.41442  0.85943  0.75000
      0.82080  0.37914  0.25000
      0.12972  0.59495  0.75000
      0.51237  0.32643  0.03597
      0.68259  0.78406  0.95699
      0.70447  0.06322  0.04214
      0.25132  0.91515  0.95880
      0.90062  0.60499  0.03995
      0.06023  0.37276  0.95959
      0.68259  0.78406  0.54301
      0.51237  0.32643  0.46403
      0.25132  0.91515  0.54120
      0.70447  0.06322  0.45786
      0.06023  0.37276  0.54041
      0.90062  0.60499  0.46005