材料科学中的数据挖掘：晶体图神经网络解读与代码解析
本研究提出了一种基于晶体图神经网络（CGNN）的模型，用于预测材料的形成能、晶胞体积、带隙和总磁化强度。通过将材料数据转换为晶体图，模型能够有效捕捉晶体结构与材料性质之间的关系。实验结果表明，该模型在OQMD数据库上的预测准确性接近基于密度泛函理论（DFT）计算的结果，且在分类材料性质（如金属、绝缘体、磁性）时表现优异。研究还展示了晶体图神经网络在材料科学中的潜在应用，为开发高效材料设计工具提供了理论依据。该模型通过预测材料性质和分类功能，为材料科学和工程领域提供了新的研究方向。
©PaperWeekly 原创 · 作者｜张玮玮
学校｜东北大学硕士
研究方向｜情绪识别
论文标题：
Crystal Graph Neural Networks for Data Mining in Materials Science
论文链接：
https://storage.googleapis.com/rimcs_cgnn/cgnn_matsci_May_27_2019.pdf
代码链接：
https://github.com/Tony-Y/cgnn
该数据库规模较大，提供 Python API 和晶体结构可视化功能，支持通过 MySQL 导入数据：OQMD [2]
晶体材料的性质通常通过计算得出，这些计算基于晶胞内的原子数目、原子分子坐标以及晶格常数，使用密度泛函理论（DFT）进行。机器学习方法在加速新材料设计方面日益普及，其在预测材料性能方面的准确性与从头算计算相媲美，但计算速度却快了数量级。最近，有人提出了一种称为晶体图的新表示方法[1]，该方法利用平衡键距作为空间信息，成功通过晶体图卷积神经网络预测了材料的体积特性。
预备知识
1.1 晶格常数（点阵常数 ）
晶胞的三个棱长及其间的夹角构成了6个参数。根据晶格常数的不同，可以将晶体划分为七大晶系。在本文中，点阵矢量被定义为。
1.2 原子坐标
原子坐标有两种表示形式：分数坐标和笛卡尔坐标。分数坐标是将点阵矢量视为单位矢量来确定原子位置，而笛卡尔坐标则是基于直角坐标系和斜角坐标系的绝对坐标体系，其计算需考虑点阵参数的大小。本文采用的是原子的分数坐标表示。笛卡尔坐标与分数坐标之间可以相互转换，在定义的晶胞中，第n个原子的笛卡尔位置计算为：(xₙ, yₙ, zₙ) = (aXₙ, bYₙ, cZₙ)，其中晶胞体积V = abc√(1 - α² - β² - γ² + 2αβ + 2αγ + 2βγ)必须大于0。
1.3 晶体性质
晶体结构决定了晶体材料的总能量。此外，DFT方法可以用来近似计算，其中，原子数目为。其近似能量表达式为：，通过最小化总能量，可以确定平衡结构。
平衡能的最小值被表示为E_{min}。DFT计算提供了多种体积特性指标，包括形成能的值、单位单元体积的平均值、带隙的大小以及总磁化强度的分布。
1.4 数据集
现有理论材料数据库涵盖了无机晶体材料的松构体结构及其性能参数。The Materials Project（MP）数据库包含86,000条目数据，其中大部分与实验观察到的结构信息相关。The Open Quantum Materials Database（OQMD v1.2）则包含563,000条目数据，这些条目数据主要涉及实验和假设化合物的相关研究。
晶体图神经网络
晶体图神经网络（CGNNs）开发了一种比例保持不变的图协调器（Crystal Graph Coordinator）来构建晶体结构，用于训练CGNN模型于OQMD数据集。该模型对每一种测试材料的形成能、单元体积、能隙和总磁性强度等体积特性进行了预测分析，其预测误差均低于数据库中的相应误差。预测的能隙和总磁性强度被用于金属绝缘体与非磁-磁材料的分类任务。
该网络架构，如图1所示，由CG-Gen模块构建，其生成的原子序数被视为自变量，能够生成平衡结构下的晶体图。该模型旨在预测晶体间的平衡特性，适用于多任务学习场景。该变体则在特定条件下推导出晶胞形状和原子分数坐标，从而实现平衡结构的构建。通过预测的晶胞体积，我们能够对晶胞进行缩放处理，以获得理想的平衡晶胞结构。
2.1 The Crystal Graph Coordinator
为了训练 CGNN 和 CG-Gen 晶体图，该框架设计了一种平衡结构。在晶体学研究中，泰森多边形被用于确定原子的最近邻居，参考文献 [1]对此方法进行了应用。CGNN 通过引入一种晶体图协调器，该协调器利用聚类技术将原子间的间距进行连接。为了确保晶体图在单位单元的任何均匀变形下保持不变性，该框架采用了体积校正因子和缩放距离。伪代码如图所示。
源代码如下所示：
def get_neighbors(geom):
elems = [Element.from_Z(z) for z in geom.atomic_numbers]#通过原子序数来得到元素种类；OQMD数据库共有89种化学元素，原子序数为1到83,89到94
radii = np.array([get_radius(e) for e in elems])#原子半径是从PyMatGen库获得
cutoff = radii[:,np.newaxis] + radii[np.newaxis, :]#将任意两个原子的半径相加，比如有4种原子，则返回一个4*4的矩阵。
vol_atom = (4 * np.pi / 3) * np.array([r**3 for r in radii]).sum()#将晶胞内的原子体积相加
factor_vol = (geom.volume / vol_atom)**(1.0/3.0)#体积校正因子通过晶胞体积和原子体积之和得来
factor = factor_vol * RADIUS_FACTOR
cutoff *= factor#得到截止半径
candidates = get_nbrs(geom.cart_coords, geom.lattice.matrix, cutoff)#这部分代码函数在下文提供，主要是将原子的笛卡尔坐标转换为分子坐标并且得到当前原子小于截止半径范围内的候选邻居原子。
neighbors = []
for j in range(len(candidates)):
dists = []
for nbr in candidates[j]:
i = nbr[0]#第一维表示原子索引
d = nbr[1]#第二维表示距离
r = nbr[2]
dists.append(d / cutoff[j,i])#距离归一到[0,1],以便进行后续的聚类
X = np.array(dists).reshape((-1, 1))
nnc_nbrs = get_nnc_loop(X)#得到最近邻居原子的位置索引
neighbors.append([candidates[j][i][0] for i in nnc_nbrs])
return neighbors def get_nbrs(crystal_xyz, crystal_lat, R_max):
A = np.transpose(crystal_lat)
B = np.linalg.inv(A)
crystal_red = np.matmul(crystal_xyz, np.transpose(B))#笛卡尔坐标系转换为分数坐标系
crystal_nbrs = pbc.get_shortest_distances(crystal_red, A, R_max,
crdn_only=True)#得到当前原子小于截止半径范围内的候选邻居原子
return crystal_nbrs得到的晶体图如下例所示：
2.2 Crystal Graph Neural Networks
作者对模型结构的介绍可以见 Architectures-CGNN（tony-y.github.io）[3]
CGNN 体系结构从原子序数 的嵌入层开始构造，得到初始隐藏状态： 。
通过T个卷积块，可以生成一系列更高维的隐藏特征：
本文开发了EdgeNet层的三个变体版本：基础版本的EdgeNet层、优化版本的EdgeNet层以及综合汇总版本的EdgeNet层。
门控卷积表示为M_{kc}，MFCNet各层采用权值矩阵。graph-level表示为是由所有隐藏层状态。在每一步中，隐藏状态采用Gated Pooling机制被聚合为：
然后，graph-level 状态表示为 被加权平均为：
将 graph-level MFCNet 的输出结果作为线性回归模型的输入，用于预测目标值。
class GGNN(nn.Module):
"""
Gated Graph Neural Networks
5. Nodes -> Embedding -> Gated Convolutions -> Graph Pooling -> Full Connections -> Linear Regression
"""
def __init__(self, n_node_feat, n_hidden_feat, n_graph_feat, n_conv, n_fc,
activation, use_batch_norm, node_activation, use_node_batch_norm,
edge_activation, use_edge_batch_norm, n_edge_net_feat, n_edge_net_layers,
edge_net_activation, use_edge_net_batch_norm, use_fast_edge_network,
fast_edge_network_type, use_aggregated_edge_network, edge_net_cardinality,
edge_net_width, use_edge_net_shortcut, n_postconv_net_layers,
postconv_net_activation, use_postconv_net_batch_norm, conv_type,
conv_bias=False, edge_net_bias=False, postconv_net_bias=False,
full_pooling=False, gated_pooling=False,
use_extension=False):
super(GGNN, self).__init__()
#设置激活函数，论文共设置了softplus,ssp,elu,relu,selu,celu六种激活函数
act_fn = get_activation(activation)
if node_activation is not None:
node_act_fn = get_activation(node_activation)
else:
node_act_fn = None
if edge_activation is not None:
edge_act_fn = get_activation(edge_activation)
else:
edge_act_fn = None
postconv_net_act_fn = get_activation(postconv_net_activation)
#EdgeNet层的三个变体
if n_edge_net_layers < 1:
edge_nets = [None for i in range(n_conv)]
else:
edge_net_act_fn = get_activation(edge_net_activation)
if use_aggregated_edge_network:
#AggregatedEdgeNetwork
edge_nets = [AggregatedEdgeNetwork(n_hidden_feat, n_edge_net_feat,
n_edge_net_layers, cardinality=edge_net_cardinality,
width=edge_net_width, activation=edge_net_act_fn,
use_batch_norm=use_edge_net_batch_norm,
bias=edge_net_bias,
use_shortcut=use_edge_net_shortcut)
for i in range(n_conv)]
elif use_fast_edge_network:
#FastEdgeNetwork
edge_nets = [FastEdgeNetwork(n_hidden_feat, n_edge_net_feat,
n_edge_net_layers, activation=edge_net_act_fn,
net_type=fast_edge_network_type,
use_batch_norm=use_edge_net_batch_norm,
bias=edge_net_bias,
use_shortcut=use_edge_net_shortcut)
for i in range(n_conv)]
else:
#Original EdgeNet Layer
edge_nets = [EdgeNetwork(n_hidden_feat, n_edge_net_feat,
n_edge_net_layers, activation=edge_net_act_fn,
use_batch_norm=use_edge_net_batch_norm,
bias=edge_net_bias,
use_shortcut=use_edge_net_shortcut)
for i in range(n_conv)]
if n_postconv_net_layers < 1:
postconv_nets = [None for i in range(n_conv)]
else:
postconv_nets = [PostconvolutionNetwork(n_hidden_feat, n_hidden_feat,
n_postconv_net_layers,
activation=postconv_net_act_fn,
use_batch_norm=use_postconv_net_batch_norm,
bias=postconv_net_bias)
for i in range(n_conv)]
self.embedding = NodeEmbedding(n_node_feat, n_hidden_feat)#节点嵌入层，使用线性层
self.convs = [GatedGraphConvolution(n_hidden_feat, n_hidden_feat,
node_activation=node_act_fn,
edge_activation=edge_act_fn,
use_node_batch_norm=use_node_batch_norm,
use_edge_batch_norm=use_edge_batch_norm,
bias=conv_bias,
conv_type=conv_type,
edge_network=edge_nets[i],
postconv_network=postconv_nets[i])
for i in range(n_conv)]#门卷积层
self.convs = nn.ModuleList(self.convs)
if full_pooling:
n_steps = n_conv
if gated_pooling:
self.pre_poolings = [GatedPooling(n_hidden_feat)
for _ in range(n_conv)]
else:
self.pre_poolings = [LinearPooling(n_hidden_feat)
for _ in range(n_conv)]
else:
n_steps = 1
self.pre_poolings = [None for _ in range(n_conv-1)]
if gated_pooling:
self.pre_poolings.append(GatedPooling(n_hidden_feat))
else:
self.pre_poolings.append(LinearPooling(n_hidden_feat))
self.pre_poolings = nn.ModuleList(self.pre_poolings)
#门pooling层
self.pooling = GraphPooling(n_hidden_feat, n_steps, activation=act_fn,
use_batch_norm=use_batch_norm)
#多层全连接神经网络
self.fcs = [FullConnection(n_hidden_feat, n_graph_feat,
activation=act_fn, use_batch_norm=use_batch_norm)]
self.fcs += [FullConnection(n_graph_feat, n_graph_feat,
activation=act_fn, use_batch_norm=use_batch_norm)
for i in range(n_fc-1)]
self.fcs = nn.ModuleList(self.fcs)
#线性回归
self.regression = LinearRegression(n_graph_feat)
if use_extension:
self.extension = Extension()
else:
self.extension = None
def forward(self, input):
#模型顺序如下：Nodes -> Embedding -> Gated Convolutions -> Graph Pooling -> Full Connections -> Linear Regression
x = self.embedding(input.nodes)
y = []
for conv, pre_pooling in zip(self.convs, self.pre_poolings):
x = conv(x, input.edge_sources, input.edge_targets)
if pre_pooling is not None:
y.append(pre_pooling(x, input.graph_indices, input.node_counts))
x = self.pooling(y)
for fc in self.fcs:
x = fc(x)
x = self.regression(x)
if self.extension is not None:
x = self.extension(x, input.node_counts)
return x
实验结果
研究使用的数据集包含从 OQMD 中收集了 561,888 个条目的相关数据，其中包括了各条目对应的松弛结构、形成能、晶胞体积、带隙和总磁化率等关键参数。该模型分别对形成能、晶胞体积、带隙这三个物理性质进行了预测任务的设置，并通过实验验证了其预测的准确性。具体预测结果如下：
测试集中，非磁体样本构成30,299个，磁体样本构成25,629个。非磁体组的总磁化率几乎为零，其占总磁化率的比重为54.2%。而磁性组的平均值为X /atom，标准偏差为Y /atom，极大值为Z /atom。与带隙数据类似，全题磁化数据在零点处呈现非负偏置特征，然而磁数据表现出显著的正偏度，其值为1.66，且峰值度为4.72，这表明数据分布具有重尾特征。基于金属绝缘体与非磁-磁材料的分类结果，
结果与讨论
对于OQMD测试集中的每一种材料，P-CGNN模型准确地预测了其形成能、单元体积、带隙和总磁化强度，并将其划分为金属和绝缘体、非磁体和磁体。值得注意的是，预测未采用空间距离键长信息。此外，CGNN模型其预测精度与基于DFT的材料性质预测相媲美。
当CG-Gen从适当材料数据库中几乎完全掌握了晶体图上的合成化学知识时，其便能够实现材料数据挖掘任务。晶体图形序列是由对应的图形转换序列生成的，而整个图形转换序列集合则构成了物质基因组，类似于人类基因组存储了核酸序列的完整信息。
参考文献
Crystal-graph-based convolutional neural network models are proposed to achieve accurate and interpretable predictions of material properties. This innovative approach has been validated through extensive computational experiments, demonstrating its effectiveness in analyzing complex material structures.
[3] https://tony-y.github.io/cgnn/architectures/
更多阅读
[
该文章通过分享这一优质内容，为读者提供了深入的见解与建议。
[
该文章通过...原式推导出...结果，为解决...问题提供了新的思路。该研究方法具有显著的效果，能够有效提升...的效率。研究者建议在实际应用中结合...理论进行深入分析，以进一步验证其可行性。
[
#投 稿 通 道#
让你的论文被更多人看到
如何能让更多的优质内容通过更短的路径到达读者群体，从而降低读者寻找优质内容的成本？答案就是：那些不被你认识的人。
某些时候，总会有一些你不认识的人，具备你所需的某些知识储备。PaperWeekly 或许能成为这样一座桥梁，促进不同学术背景和研究方向的学者之间的交流，激发更多的创新火花。
PaperWeekly 欢迎高校实验室或个人，在我们的平台上分享各类优质内容。这些内容可以是最新论文解读、学习心得以及技术干货。我们的宗旨很简单，就是让知识真正流动起来。
???? 来稿标准：
稿件为个人原创作品，请作者在来稿中明确标注个人信息（包括姓名、所属学校或工作单位、学历或职位以及研究方向）。
• 如果文章并非首发，请在投稿时提醒并附上所有已发布链接
• PaperWeekly 默认每篇文章都是首发，均会添加“原创”标志
???? 投稿邮箱：
• 投稿邮箱：hr@paperweekly.site
• 所有文章配图，请单独在附件中发送
• 请留下即时联系方式（微信或手机），以便我们在编辑发布时和作者沟通
????
现在，在**「知乎」** 也能找到我们了
进入知乎首页搜索**「PaperWeekly」**
点击**「关注」** 订阅我们的专栏吧
关于PaperWeekly
PaperWeekly 是一个致力于推荐、深入分析、分享观点以及报道人工智能前沿论文成果的学术平台。如果你是研究人员或AI从业者，欢迎在公众号后台点击「交流群」，通过小助手，你可以加入 PaperWeekly 的专业交流圈子。
