Complete and Efficient Graph Transformers for Crystal Material Property Prediction

Thank you for your recognition of the presentation and the problem we focused on.

1. This method is not universal as it cannot handle crystals with lattice containing (0,0,0). The paper mentions that “Third, we select the duplicate with the next smallest distance ||e_ii3||_2, and verify that eii3, eii1, and eii2 are not in the same plane in 3D space, and repeat until eii3 is found”, but if lattice containing (0,0,0), then eii3 cannot be found. How to deal with the crystal whose lattice contains (0,0,0)? Do the datasets in experiments contain this type of data? In theory, the method proposed by the author cannot handle this situation. Why was it not explained in experiments?

Thank you for this question. For a 3D crystal structure, following the traditional and widely-used definition, if we have a third lattice vector $l_3$ = (0, 0, 0), it means for a given atom i in the cell with position p_i, there will have infinite number of atoms at position $p_i$, due to $p_i' = p_i + k_3 \* l_3 = p_i$. Thus, this scenario will not happen in reality. If you are talking about purely 2D crystals with a single layer and strictly all z=0, then only the first two lattice vectors are needed to achieve completeness.

For 2D crystals not strictly have all z=0, the traditional way to represent the structure is setting $||l_3||_2 > 10$ Angstroms and not in the same plane of $l_1$, $l_2$, and the proposed method can naturally handle this case.

It is worth noting that our experiments include 2D crystals - matbench jdft2d as shown in the main paper. To make this clear, we revised the paper accordingly and added more details about jdft2d in the experimental section and the appendix.

2. The experiment is not comprehensive enough. The most significant contribution of this paper is that it proposes a complete representation, but the final experiment only verifies that it can better predict without verifying the completeness of the representation. If possible, it is recommended that the author add some experiments that demonstrate the completeness of the representation.

Thank you for this question. In the original version, we put the completeness verification experiments in Appendix A.4. We do believe in the importance of completeness of our proposed representations, so we follow your suggestion and further conduct completeness verification of our graph representations on a more challenging T2 crystal dataset in which the largest crystal cell contains more than 700 atoms. It can be seen that the proposed graph representations can fully reconstruct the super large crystal structures with only k=25 (where k means k nearest neighbors). Specifically, We convert crystal structures to the proposed SE(3) invariant and SO(3) equivariant crystal graphs and then use reconstruction steps shown in Appendix A.4 to recover the crystal structures from the graphs. We then compare the structure differences between the recovered ones and the input ones using Pymatgen. As shown in the below table, the input crystal structures are fully recovered, with limited structure error (less than $10^{-6}$).

According to your suggestion, we further moved the completeness verification experiments from the Appendix to the main paper, to enhance the clarity.

[1] T2 dataset: Pulido, Angeles, et al. "Functional materials discovery using energy–structure–function maps." Nature 543.7647 (2017): 657-664.

3. What is the role of SE(3) invariant and SO(3) equivariant? In what scenarios should eComFormer be used, and in what scenarios should iComFormer be used? What are the reasons for iComFormer outperforming eComFormer in most cases?

For the tasks of scalar property prediction, SE(3) and SO(3) can both be used. However, for higher tensor order crystal property predictions, only SO(3) will be suitable. We suggest the usage of iComFormer for scalar property predictions and the usage of eComformer for the development of higher tensor order predictions, which is also one of our future focus.

The reason that iComFormer performed a little better than eComFormer (actually only marginal in many cases, because both of them use complete graphs) is due to the equivariant tensor product layers may require more training samples and are harder to train than invariant features because equivariant features cannot support nonlinear mappings.

We believe your suggestions enhanced the clarity of this paper.

Yours sincerely, authors

Thank you for the author's response, and most of my doubts and concerns have been addressed.

Here are some of my new questions.

The role of eComformer in predicting the properties of tensor ordered crystals has not been reflected in this work, so it lacks the persuasiveness to consider this as the difference between the two methods. I would like to know the role of eComformer in this work and the differences between these two methods in the tasks involved in this paper.

In addition, I have another question that I would like to discuss with the author. Given that I am not an expert in the field of materials, I am curious to know whether graphene, being a two-dimensional structure in the real world, encounters the issue of the lattice containing (0,0,0).

4. Is it possible to provide some qualitative experimental results (visualization) of SE(3) invariant property and/or SO(3) equivariant property, to prove the graph representation and the proposed Transformer network indeed learns such properties? It is not mandatory, but just out of curiosity.

Sure, we provided further experimental verification of the SE(3) invariance of the proposed eComformer and iComformer for scalar properties. We use the MP dataset test set, and calculate the prediction differences of every test sample when the structure is rotated by R = [[0, -1, 0], [1, 0, 0], [0, 0, 1]], and found the differences before and after this rotation for all test samples are very small (less than 10^{-5} introduced by numerical errors). We further show that when a crystal structure ZnS (mp-560588) without inversion symmetry is reflected by R = [[-1, 0, 0], [0, -1, 0], [0, 0, -1]], the constructed graphs are totally different. For iComformer, the three lattice angles (cosine values) before inversion between atom 0 and atom 3 is (-0.8527, -0.8527, -0.1751) and is (0.8527, 0.8527, 0.1751) after inversion. For eComformer, the edge features $e_{ij}$ for atoms 0 and 3 before inversion is (3.8058, -2.1973, 0.7815) and changes to (-3.8058, 2.1973, -0.7815) after inversion, while lattice representation remain the same [-3.8058,0.,0.], [-1.9029,3.2959,0.], [ 0., 0., -6.2380], resulting in different SO(3) graphs.

Thus, iComformer and eComformer are SE(3) invariant when predicting scalars.

5. In Table 5, the size of eComFormer is larger than iComFormer by two to three times. Any explanation or reason for this model size?

The parameter numbers of eComformer mainly come from tensor product layers that are more expensive than normal message passings. We further conducted a stress test by decreasing the number of channels in equivariant tensor product layers by half, and the eComformer will only have 5.6 M parameters and is close to iComformer (4.1 M), the performance of the reduced eComformer is 29.1 meV/atom, which is very close to the full model (28.4 meV/atom), and still the SOTA performance beyond all previous works.

6. The paper spends a lot of space on background introduction and graph representation description, leading to just 1.5 pages for Transformer network design. It lacks the description of the rationale behind the network design, especially for SO(3) message passing, which causes a bit confusion in finding enough novelty of the proposed network.

Thank you for this suggestion, as discussed in the answer to question 3, The aim of the proposed eComformer message passing will be, how we can make the best use of these $e_{ij}$ to capture rich geometric information without breaking SO(3) equivariance during message passing.

Following your suggestion, we revised SO(3) message passing section and Appendix A.5 to emphasize these points and the design considerations. We believe your suggestions enhanced the clarity of this paper.

Yours sincerely, authors

Thank you for your recognition in terms of motivation, sound graph construction, and experimental results. We provide point-to-point responses for your mentioned questions as below.

1. Is it just the way of lattice representation construction (Section 3.2 and Section 3.3) that introduces SE(3) invariant and SO(3) equivariant properties, or the proposed Transformer network (Section 4) have such properties? In other words, if the proposed Transformer network is applied to an ordinary graph (say a k-NN graph), do the SE(3) invariant and SO(3) equivariant properties still hold?

Thank you for your question. The constructed graphs introduce SE(3) invariant and SO(3) equivariant properties, while the proposed Transformer networks maintain such properties. The iComformer and eComformer are specially designed for corresponding SE(3) invariant and SO(3) equivariant graphs, and cannot take ordinary graphs as input (ordinary graphs do not have that 3 lattice angles). The whole pipeline satisfies SE(3) invariance and they need to be used together.

2. It seems that the node-wise transformer layer in SE(3) invariant message passing (Section 4.1) has no significant difference from the conventional graph Transformer (or graph attention network) [1]. Though the detailed design may be different, the general structure/scheme is similar. Correct me if I'm wrong.

We respectfully disagree with this. The major differences are as follows.

First, GAT only considers node features and can not model 3D objects. The attention mechanism is $\alpha_{ij} = W(h_i, h_j)$, where $\alpha_{ij}$ is the attention weight for node pair i and j, and $h_i$ is the node feature for i.

Second, according to $\alpha_{ij} = W(h_i, h_j)$, GAT has node-wise attention and can only model ordinary graphs with a single edge between nodes, while crystal graphs are multi-edge graphs.

Different from GAT, SE(3) invariant message passing captures 3D information including Euclidean distances and lattice angles through edge-wise attention and is specifically designed for crystals.

3. It is unclear why edge feature needs to be embedded using spherical harmonics, in second paragraph of Section 4.2. A justification is needed. What is tensor product layer and why need it here? Could you provide a concise and intuitive explanation?

We merge these two questions together since they are closely related. Let's begin with the features in SO(3) crystal graph and then discuss why TP layers are needed and why embed edge feature vectors into spherical harmonics. The edge features in SO(3) graph are $e_{ij} = p_j - p_i \in R^3$, which will change accordingly when you rotate the coordinate system by a rotation matrix. The aim of the proposed eComformer message passing will be, how we can make the best use of these $e_{ij}$ to capture rich geometric information without breaking SO(3) equivariance.

To keep the equivariance property and finally obtain invariance predictions, there are limited operations that can be done to $e_{ij} = p_j - p_i \in R^3$. For example, one cannot use a nonlinear function like relu to $e_{ij}$. The limited operations one can use are norm(), scale by a scaler, or tensor product. However, the first two, norm() or scale, cannot capture angular information and only contain Euclidean distance information, while TP layers have the potential to capture this if properly designed like our proposed eComformer.

Concretely, the Tensor product layer takes two inputs with different rotation orders, where a scalar belongs to rotation order 0, a 3D vector belongs to rotation order 1, and rotation order 2 vectors are of length 5. For example, when you use the tensor product with two 3D vectors, the operations have been done in TP($v_1$, $v_2$) will be, (1) calculate the tensor product $v_1 \times v_2$ with shape 3*3, and then decompose this 3*3 matrix into one scalar, one vector, and one 5D vector, through Clebsch–Gordan coefficients. The resulting scalar, vector, and high rotation order vector will rotate accordingly through Wigner_D conditioned on the rotation of the whole system (e.g. Wigner_D(R) for l=1 3D vectors will be exactly R). And more mathematical details can be found in [1].

The reason to convert $e_{ij}$ to spherical harmonics is to feed it properly to TP layers, by normalizing the vector and converting it to one 3D vector with rotation order 1 (the direction of the vector) combined with a scalar with rotation order 0 (the length of the vector).

Thus, we use Tensor product layers to process SO(3) equivariant crystal graphs, and prove that the design of eComformer considers all two-hop bond angle information in Appendix A.5. We also revised the paper accordingly in Sec. 4.2 to enhance the clarity.

[1] Geiger, Mario, and Tess Smidt. "e3nn: Euclidean neural networks." arXiv preprint arXiv:2207.09453 (2022).

Dear Reviewer ZKHH,

As the rebuttal period is approaching its end, we would be grateful if you could take a moment to review our responses. Your further insights will be immensely beneficial for us to improve our work and for the committee to make a well-informed decision.

We understand that reviewing is a demanding task and truly appreciate the time and effort you have already invested in reviewing our work. In response to your valuable feedback, we have provided a detailed rebuttal addressing the questions and concerns raised.

Please feel free to reach out if there are any additional clarifications or information you might need from our end.

Thank you once again for your valuable time and consideration. We look forward to your response.

Best regards, authors

Dear reviewer ZKHH,

Since the deadline for the author-reviewer discussion is approaching soon, we would greatly appreciate it if you could provide your feedback within the remaining hours of the discussion period.

Your response will be instrumental in helping us refine our submission and contribute more effectively to the field.

We understand that you have a demanding schedule, and we appreciate the time and effort you dedicate to the review process. If there are any specific areas of our submission that you feel require further clarification or improvement, please do not hesitate to share your thoughts. We are fully committed to making the necessary revisions and providing additional information as needed.

Yours sincerely, authors

Continued answer for question 15.

-> Additionally, these descriptors cannot deal with dynamic crystal structures. For PDD (2,4), when two atom switches positions, the crystal changes but PDD remains the same, and this case will occur in dynamic systems during trajectories. For example, graphene and boron nitride (two famous 2D materials), both have hexagonal lattice, with close lattice constant of 2.46 A for graphene, and lattice constant of 2.54A. If the lattice constants are changed to be the same under hypothetical dynamic states, then these two completely different crystals will have exactly the same PDD representations, and cannot be distinguished. Our proposed method can still distinguish this case, as we do not rely on any assumption about the stability of the crystal while proving completeness.

-> Furthermore, these Euclidean distances-based descriptors cannot distinguish chiral crystals with different crystal properties.

[10] Gong, Sheng, et al. "Examining graph neural networks for crystal structures: limitations and opportunities for capturing periodicity." Science Advances (2023)

Yours sincerely, authors

Below, we provide a brief overview of each work and its traditional applications, followed by an explanation of how our research diverges in its focus and objectives.

• Niggli (1927): Paul Niggli’s contributions in 1927 laid the groundwork for modern crystallography through the development of the Niggli reduction algorithm. This algorithm standardizes crystal descriptions, making it easier to compare different crystal structures. Traditionally, Niggli's work has been pivotal in crystallographic data standardization, aiding in the systematic classification and comparison of crystalline materials.
• Lawton (1965): Lawton's work in 1965 further extended the understanding of crystal symmetry and its practical applications. Lawton's contributions are particularly noted in the realms of inorganic compound charaterization, where understanding symmetry plays a crucial role in deciphering the physical and chemical properties of materials.
• Parthe (1987): Parthe’s work in the late 1980s focused on the systematic description of crystal structures, contributing significantly to the field of structural crystallography. His methods have been extensively applied in the identification and categorization of complex inorganic structures. Parthe developed STRUCTURE TIDY code to help standardize crystal structure data and make it easier to identify identical or nearly identical structures by two standarization parameters \Gamma and CG. Optionally, this code can transform structures into Niggli's reduced cells and make it easy for structural comparison and finding out symmetry elements.

While we recognize the importance and broad applications of these classical methods in crystallography, our research primarily focuses on the predictive modeling of material properties using advanced graph neural networks. These classical methods aimed at structural standardization and comparison. However, there is a huge gap from standardized (A, P, L) to crystal property prediction, which is exactly DFT algorithms and our proposed method want to resolve.

In conclusion, while our research is informed by the rich history of crystallography, including the foundational work by Niggli, Lawton, and Parthe, it diverges in its application of computational methods for crystal property prediction. We believe that acknowledging these distinctions offers clarity on how our work contributes uniquely to the field.

We hope this response adequately addresses your query and demonstrates the contextual understanding of our research within the broader landscape of crystallography.

15. What results in the paper are stronger than the past work below?

(1) continuity under perturbations in Theorem 4.2 and generic completeness of the density fingerprint in Theorem 5.1 from Edelsbrunner et al (SoCG 2021), (2) continuity under perturbations in Theorem 4.3 and generic completeness with an explicit reconstruction in Theorem 4.4 from Widdowson et al (NeurIPS 2022), (3) full completeness in Theorem 10 from Anosova et al (DGMM 2021), extended in arxiv:2205.15298, (4) 200+ billion pairwise comparisons of all periodic crystals in the Cambridge Structural Database completed within two days on a desktop, see the conclusions in Widdowson et al.

We respectfully disagree with this. The (1,2,3,4) you mentioned are designed to identify the similarity of stable crystal structures with vector or matrix-based descriptors. Our target is to predict the properties of crystal structures, which can be stable crystals or dynamic snapshots. The clear advantages of this paper beyond (1,2,3,4) are listed below.

-> 1. These complete descriptors can only be applied to stable crystal structures and do not consider atom types. This makes them not suitable to do property predictions. For the property prediction task, we directly compared with [2] based on your mentioned (2,4), and found that when (2,4) is used to predict properties, it is no longer complete and uses graphs similar to CGCNN. While through extensive experiments, we have shown significant performance gaps beyond CGCNN, and through direct comparison with [2], we show huge improvements (more than 40% decrease in MAE).

-> 2. Also, [2] based on (2, 4) cannot achieve completeness due to practical issues when doing property prediction. To achieve completeness of (2,4), the k, where k means k-nearest neighbors, needs to be predefined to be larger than the maximum number of atoms in the test set. This is against the logic of looking into the test set or assuming the largest cell the method will encounter. For the T2 dataset[8], k needs to be larger than 700, which results in 700 * n edges to achieve completeness. However, as we show in the revised paper, our method can fully reconstruct crystal structures in the T2 test set by only using k=25.

continued in Part 6.

Continued answer for question 11: The advanced electronic structure methods such as DFT become more and more reliable for property predictions with evaluated error shown in the JARVIS paper [9]. However, their computational cost is high, which motivates many AI/ML works in the field of force field, including our proposed method towards accurate and fast prediction of materials properties and providing a short list of computationally screened potential candidates for experimentalists, thereby accelerating the discovery of new materials with desired properties.

Last but not least, the method proposed here can be applied to both the computational dataset and experimental dataset if available.

12. If initial property computations were too slow and most likely iterative, why not stop their computations earlier instead of writing another program simulating the simulations?

If the reviewer's question is about the curation of the dataset itself, the iterative calculations (more precisely self-consistent field calculations) are automatically stopped after the desired convergence is reached (e.g. when total energy difference or maximum residual atomic force is less than a given small threshold). Such iterative calculations need to be carefully done, otherwise, the calculated property will have large intrinsic errors. In other words, the initial calculations, though time-intensive, need to be carefully performed despite being slow, which is essential for generating accurate and high-fidelity data. These datasets such as MP and Jarvis serve as a robust foundation for training ML models, ensuring their predictive power and reliability (see the quality and accuracy addressed in the response to your question 11).

If the reviewer's question is about materials modeling and simulations, it is worth noting that the physics, chemistry, and materials communities face concurrent challenges. But that's exactly why AI/ML approaches can be leveraged to advance the field of natural science. To be more specific, the materials world is complex and multiscale in nature [Handbook of Materials Modeling, S. Yip, Springer] which can span a huge length scale (10 orders of magnitude from Angstroms to meters) and a vast time scale (18 orders of magnitude from attoseconds to seconds/hours). DFT-based electronic structure calculations scale as $N^3$ where $N$ is the number of electrons/atoms for each electronic self-consistent field / iterative loop. So, for a single material/molecule, the computational cost of DFT calculations is reasonable and manageable, which can be typically done within a day or a few days, depending on (1) specific properties of interest (total energy/atomic forces, quasiparticle bandgap, exciton energy), (2) configuration setup (hybrid vs local/semilocal exchange-correlation energy functionals, whether or not considering spin-orbit coupling, etc.), and (3) computational resources and codes (CPU vs CPU+GPU parallelization, plane-wave basis set vs localized/polarization basis set vs dual basis set). However, it becomes very slow for materials/molecules with >1,000 atoms accompanied by the rapidly increased memory consumption. It becomes even worse for ab initio molecular dynamics calculations where very often thousands or more ionic steps are required for statistically meaningful physical properties, despite that such ab initio simulations are highly valuable for understanding and predicting many physical and chemical processes such as chemical reactions and catalyst design. By curating high-quality and consistent datasets of small systems and developing ML models such as iComformer and eComformer here, it will be possible to bridge both length-scale and time-scale bottlenecks mentioned in the beginning and advance the scientific understanding and predictions that are not possible with conventional approaches, as demonstrated in many recent works (e.g. Physical Review Letters 98, 146401 (2007); Physical Review Letters 104, 136403 (2010); Advanced Materials 31, 1902765 (2019); Nature Communications 11, 2654 (2020); Chemical Reviews 121, 10142 (2021)).

13. How many hidden parameters and CPU hours were used for producing the experimental results?

As shown in the main paper, we use 4M params for iComformer. The training process is conducted on a single RTX A6000 GPU, where iComformer with 3 node layers takes 9.6 hours to train and eComformer 16 hours (much faster than previous SOTA ALIGNN).

14. Did the authors know about the classical results in crystallography cited above, starting from Niggli (1927), Lawton (1965), and Parthe (1987)?

Thank you for your insightful query regarding our knowledge and consideration of seminal works in crystallography, specifically those by Niggli (1927), Lawton (1965), and Parthe (1987). We acknowledge the foundational importance of these studies in the field of crystallography and their various applications over the decades. (continued in Part 5).

9. Since the paper title promised complete transformers, have the authors found any geometric duplicates in the used datasets?

The application of filtering duplicate crystals is not our focus. However, we do the following transformations to convert our proposed graphs to matrix form (Comformer-matrix) and do the screening. Following [1], we convert the crystal graph into a $n \times k \times 4$ matrix, where n is the number of atoms in a cell, k is the k-nearest neighbors, and 4 means we have 4 invariant features for every neighbor, including Euclidean distances, and 3 lattice angles. After this transformation, we can compute a $k*4$ invariant vector for this structure, and then use MAE between vectors to determine the structural differences between crystals (following AMD). By doing so, we screened the Jarvis dataset, and identified more than 100 duplicate pairs including ('JVASP-86097', 'JVASP-86697'), ('JVASP-21210', 'JVASP-1047'), ('JVASP-86322', 'JVASP-87063'), using threshold $10^{-6}$. The whole process, including graph construction, matrix transformation, and screening, takes less than 15 minutes.

Additionally, contradicting to the reviewer’s claim, the two materials mentioned by the reviewer (mp-568619 and mp-568656), both of which are silicon carbide (SiC), are NOT duplicates. Rather, they are completely DIFFERENT, due to their different stacking faults or different sequences along the z direction, known as SiC polymorphs (https://en.wikipedia.org/wiki/Polymorphs_of_silicon_carbide). This commonly happens in SiC and many other close-packed metals and semiconductor crystals (e.g. titanium Ti, gallium arsenide GaAs, etc.), which are well-known in the materials science and semiconductor communities. Specifically, we can use the conventional notation of A, B, and C sites to denote three different projection locations of each vertical Si-C pair on the x-y plane (https://en.wikipedia.org/wiki/Close-packing_of_equal_spheres). “A” stands for the site (0,0), B for the site (1/3, 1/3), and C for the site (2/3, 2/3) in the fractional coordinates on the x-y plane. Then, the stacking sequence for “mp-568619" is “CBACA” along the z direction, while the stacking sequence for “mp-568656" is “CBACABACAB”. Although the first 5 stacking layers in the two structures are the same (CBACA), the rest 5 SiC stacking layers in “mp-568656" is “BACAB”, different from “CBACA”, which means one cannot by simply doubling the crystal structure of mp-568619 along z direction to obtain the exactly same structure as mp-568656. In fact, there are 3 “C”, 3 “B”, and 4 “A” stacking layers in mp-568656. Thus, these two structures are completely different from each other. Another clear evidence is their different electronic band gap: mp-568619 has a band gap of 1.75 eV, while the other one has a band gap of 2.01 eV, which otherwise, based on group theory, would have been the same if these two structures were simply related by doubling the structure. Additionally, the MAE error between two crystal structures calculated by ComFormer-matrix is 0.093 (>> $10^{-6}$), indicating they can be distinguished by our ComFormer model.

It is worth noting that we are focusing on crystal property prediction, but not this task.

10. The dataset details in section A.5 could include the most important information: whether crystals are simulated or experimental, and if the "ground-truth" properties are obtained from computer modeling or physical experiments.

Thank you for your suggestion, we added the information in the appendix A.5. 30086 crystals in MP and 18865 crystals in Jarvis are experimentally validated (observed). All crystal properties are DFT-calculated.

11. It seems that all used crystals are hypothetical with simulated properties. Is the main contribution a computer program that predicts the outputs of other computer programs?

No, first, 30086 crystals in MP and 18865 crystals in Jarvis are experimentally validated (observed).

Second, it is worth noting that in [2] and [7] suggested by the reviewer, only the MP dataset and the simulated lattice energy by force field calculation from another work [8] were used in their tests.

Third, we disagree with "a computer program that predicts the outputs of other computer programs". The advances in the development of advanced electronic structure theory and exchange-correlation energy functional have significantly improved the accuracy of computational results, approaching the chemical accuracy of ~1 kcal/mol limited by experiments (including materials sample quality, inherent resolution of experimental techniques, and errors in measurements, etc.). The advanced electronic structure methods such as DFT become more and more reliable for property predictions with evaluated error shown in the JARVIS paper [9].

[9]Choudhary, Kamal, et al. "The joint automated repository for various integrated simulations (JARVIS) for data-driven materials design." npj computational materials 6.1 (2020): 173.

4. Any graph representation of a crystal or a molecule is discontinuous because all chemical bonds are only abstract representations of inter-atomic interactions and depend on numerous thresholds on distances and angles, while atomic nuclei are real physical objects.

We respectfully disagree with this. In the graph, we have not used chemical bonds but use Euclidean distances and 2-hop angles that are continuous.

5. The invariance under permutations of atoms is missed because the paper has no words "permute" or "permutation". The permutation invariance seems lost when using angles, which depend on point ordering. Did the paper consider the invariance of crystal descriptors under permutations of atoms?

We respectfully disagree with this. To verify the permutation invariance of our proposed method, we randomly permute atom indexes in crystal structures from the Materials Project test set, and the predictions remain unchanged (the maximum property prediction difference before and after permutation is $2.24\*10^{-7}$ for Comformer without training, due to numerical error). This is because the graph construction process and message passing have not used any ordering information of atoms.

6. The property prediction makes sense only for really validated data, not for random values because neural networks with millions of parameters can predict (overfit) any data.

We show that 30086 crystals in MP and 18865 crystals in Jarvis are experimentally validated (observed). We also further conducted an experiment on 30086 crystals (pure experimentally validated dataset) and compared with [2] and [7] you mentioned as shown above in question 3. For "because neural networks with millions of parameters can predict (overfit) any data", we kindly disagree, especially for the potential crystal space which is infinite. As far as we know, there is no evidence suggesting several million parameters can predict (fit) infinite crystal space.

We assume you are talking about overfitting on the training set, but the robustness and prediction accuracy are verified on test sets of benchmarks, which are not seen by the trained models.

7. The words "crystal graph" appear a few times before section 2.2 with references but without details. In Definition 1, could the authors please explain the meaning of "a crystal graph" and the symbols "=" and "->"?

We provided the detailed crystal graph construction process with figures in Appendix 1.1, and show the detailed crystal graph construction in Sec. 3.2 and 3.3. A crystal graph in our approach consists of node features $a_i$ (atom type) for every node, and edge features including Euclidean distances between node pairs and three lattice angles for the invariant graph, and a Euclidean vector for the equivariant graph. "=" means two crystal graphs have the same node features and corresponding edge features. "A->B" means if we have A, then B is true.

8. Does a manually chosen cutoff radius (line 3 from the bottom of page 4) make any invariant incomplete because perturbing atoms can scale a unit cell larger than any fixed radius?

No. The proof of completeness is only based on the assumption that the constructed crystal graph is strongly connected. Actually, no method can achieve completeness with a k (k nearest neighbors) that will make the crystal graph not strongly connected as shown in Appendix A.8.1, including PDD[1].

Further, to show the robustness of our proposed crystal graphs in achieving completeness, we use k=25 and show that it is enough to fully reconstruct crystal structures with more than 700 atoms in a cell (T2 dataset[8] used in [2] and [7]) as shown in the revised paper in Sec. 5.1.

[1] Widdowson et al (NeurIPS 2022) -> Resolving the data ambiguity for periodic crystals. Daniel E Widdowson and Vitaliy Kurlin.

[2] Balasingham et al (arxiv:2212.11246) -> Compact Graph Representation of crystal structures using Point-wise Distance Distributions. Jonathan Balasingham, Viktor Zamaraev, and Vitaliy Kurlin.

[3] Recognizing Rigid Patterns of Unlabeled Point Clouds by Complete and Continuous Isometry Invariants with no False Negatives and no False Positives. Daniel E Widdowson and Vitaliy A Kurlin.

[4] Widdowson et al (MATCH 2022) -> Daniel Widdowson, Marco M. Mosca, Angeles Pulido, Andrew I. Cooper, Vitaliy Kurlin.

[5] Edelsbrunner et al (SoCG 2021) -> The Density Fingerprint of a Periodic Point Set. Herbert Edelsbrunner, Teresa Heiss, Vitaliy Kurlin, Philip Smith, and Mathijs Wintraecken.

[6] Anosova et al (DGMM 2021) -> Algorithms for continuous metrics on periodic crystals. Olga Anosova and Vitaliy Kurlin.

[7] Fast predictions of lattice energies by continuous isometry invariants of crystal structures. J Ropers, MM Mosca, O Anosova, Vitaliy Kurlin, AI Cooper

[8] T2 dataset: Pulido, Angeles, et al. "Functional materials discovery using energy–structure–function maps." Nature 543.7647 (2017): 657-664.

Thank you for your valuable time and comments, we provide point-to-point responses to your questions as follows.

1. The word "problem" appears only once (in the abstract), though a rigorous and explicit problem statement might have helped to understand the unresolved challenges.

Thank you for your suggestion, we revised the introduction, related works, and appendix to include the discussions with [1-7], we believe your suggestion indeed makes the paper clearer for experts in crystallography. For the detailed problem definition and challenges, we provided them in Sec. 2.1 in the original paper.

2. If we need a complete and invariant description of a periodic crystal, crystallographers solved this problem nearly 100 years ago by using Niggli's reduced cell of a lattice and then recording all atoms in so-called standard settings.

We kindly disagree with this. The discussion with Niggli's reduced cell of a lattice is included in the introduction. Niggli's reduced cell cannot be directly used to predict crystal properties, it just transforms a crystal structure to its standard form which still consists of coordinates, lattice vectors, and atom types.

3. However, all these standardizations have become obsolete in the new world of big and noisy data because the underlying lattice (not even a unit cell) of any periodic crystal is discontinuous under almost any perturbation. Moreover, atoms in any material always vibrate above absolute zero temperature, so their positions continuously change. As a result, any crystal structure with fixed atomic coordinates in a database is only a single snapshot of a potentially dynamic object. Hence the new essential requirement for any (better than the past) representations of crystals is a proven continuity under perturbations of atoms. This problem has been solved by Widdowson [1] with theoretical guarantees and practical demonstrations on the world's largest collection of materials: the Cambridge Structural Database. The underlying invariants have a near-linear time in the motif size and were used for property predictions by Ropers [7] and Balasingham [2]. Despite the authors citing one of the papers above, Definition 1 didn't follow the definitions from the past that modeled a crystal as a periodic point set not as a graph.

We agree with this point. For crystal structures with changed positions, the property will change accordingly but not drastically.

Hence, it is worth noting that the features used in our proposed graphs including Euclidean distances and two-hop angles, both are continuous under small perturbations. For the lattice, it is of great importance in determining various crystal properties as shown in a recent work [10]. All Euclidean distances, including lattice lengths, are encoded to 1/r. So for the lattice lengths a, b, c, which are usually large, changed to 2a, 4a, or 8a under perturbations, the maximum difference introduced in the representation after averaging is no larger than 1/ka (for k = 25 and a>4, this will be smaller than 0.01 which is of similar scale of a reasonable perturbation), and the proposed graph representations can capture these changes.

Moreover, we directly compared with the methods [2] and [7] you mentioned, and we found that the work [2] based on [1] you mentioned actually uses graphs to represent crystals and have very close performances with CGCNN and cannot achieve completeness due to the usage of $k=12$ (to achieve completeness, $k>200$ is needed for the materials project according to [1]) while the performance of [7] is worse than CGCNN. We conduct this comparison mainly following [2] using experimentally verified crystals from the Materials Project-2018, with 30084 crystals. We also provided more discussions with [2] and [7] in Appendix 1.3.

The performance gaps between iComformer and you mentioned [1, 2, 7] are super clear even with fewer training samples ([2] did not release the code and dataset), with more than 40% decrease in MAE. AMD [7] using the Gaussian Process needs more than 10 hours of training for 24067 samples and cannot scale to large datasets.

Thus, your mentioned works [2] and [7] do not achieve completeness when used as input to predict properties and have similar or worse performances with CGCNN, due to practical issues. We included this comparison in the appendix and will release the dataset to make it easier for the ML community to follow the insightful works you mentioned [1,2,7].

Thank you for your recognition of our method in terms of tackling practical problems in material analysis, theoretical proofs, and extensive experiments.

1. The authors mentioned learning higher-order properties in their own discussion of limitations. This could possibly be solved by the SO(3)-representation with a spherical harmonic basis.

Yes, we indeed plan to explore how to incorporate tensor field networks or TP layers to predict higher-order properties in crystals. These higher-order properties are of matrix form, for example, the dielectric tensors are rank two tensors and with shape 3*3, while elastic tensors are rank 4 tensors with shape 3*3*3*3 and are different from higher rotation order tensors generated from TFN. We plan to explore how to tackle these matrix form predictions in future works.

Yours sincerely, authors