2021年第十一届MathorCup高校数学建模
B题 三维团簇的能量预测
原题再现
原题再现
团簇，也称超细小簇，属纳米材料的尺度概念。团簇是由几个乃至上千个原子、分子或离子通过物理或化学结合力组成的相对稳定的微观或亚微观聚集体，其物理和化学性质随所含的原子数目而变化。
团簇是材料尺度纳米材料的一个概念。团簇的空间尺度是几埃至几百埃的范围，用无机分子来描述显得太小，用小块固体描述又显得太大，许多性质既不同于单个原子分子，又不同于固体和液体，也不能用两者性质的简单线性外延或内插得到。因此，人们把团簇看成是介于原子、分子与宏观固体物质之间的物质结构的新层次。团簇科学是凝聚态物理领域中非常重要的研究方向。
团簇可以分为金属团簇和非金属团簇，由于金属团簇具有良好的催化性能，因此备受关注。但由于团簇的势能面过于复杂，同时有时候还需要考虑相对论效应等，所以搜索团簇的全局最优结构（即能量最低）显得尤为困难。其中，传统的理论计算方法需要数值迭代求解薛定谔方程，并且随原子数增加，高精度的理论计算时间呈现指数增长，非常耗时。因此，目前需要对这种方法加以改进，例如：考虑全局优化算法，结合机器学习等方法，训练团簇结构和能量的关系，从而预测新型团簇的全局最优结构，有利于发现新型团簇材料的结构和性能。
请建立三维团簇能量预测的数学模型，并使用附件中的坐标和能量数据，解决下列问题。
备注：附件中数据集格式为 xyz，第一行是原子数，第二行是能量，后面是原子的三维坐标。可用文本阅读器打开，并用 VMD 等软件进行可视化。
VMD软件安装及
问题 1：针对金属团簇，附件给出了 1000 个金团簇 Au20的结构，请你们建立金团簇能量预测的数学模型，并预测金团簇 Au20 的全局最优结构，描述形状；
问题 2：在问题 1 的基础上，请你们设计算法，产生金团簇不同结构的异构体，自动搜索和预测金团簇 Au32的全局最优结构，并描述其几何形状，分析稳定性；
问题 3：针对非金属团簇，附件给出了 3751 个硼团簇 B45-的结构，请你们建立硼团簇能量预测的数学模型，并预测硼团簇 B45-的全局最优结构，描述形状；
问题 4：在问题 3 的基础上，请你们设计算法，产生硼团簇不同结构的异构体，自动搜索和预测硼团簇 B40-的全局最优结构，并描述其几何形状，分析稳定性。
整体求解过程概述(摘要)
团簇是由多个分子或原子聚集在一起的微观结构，研究团簇的全局最优结构（即能量最低）对于发现新型材料的结构和性能具有重要意义。传统的理论研究方法存在计算时间长、计算效果差等问题，而机器学习作为一个极具前景的多学科交叉领域，能够有效提高模型学习与计算的效率。因此本文针对三维团簇的能量预测和结构寻优问题，采用了多种机器学习方法进行研究。
针对问题 1，我们首先分析了传统金属团簇能量模型，将不同模型得到的能量值与样本能量值进行相关性分析，构建了含有变异因子的Sutton Chen − 能量预测模型。在此基础上改进传统 PSO 算法，引入自适应惯性权重，能够自适应调整粒子群的全局与局部搜索能力，解决了传统算法收敛速度慢、容易陷入局部最优的问题。结果表明，金团簇Au20的全局最低能量值为-1557.8251，全局最优结构为具有对称性的正四面体。
针对问题 2，我们构建了多条件约束的金团簇Au32异构体生成方法，给出了产生团簇异构体的原子结构对称性和原子间距约束条件，结合 AWPSO 算法，对异构体的绝对位置信息进行相对距离化与向量化处理，并完成了全局最优能量的预测及结构的三维重构。结果表明，金团簇Au32全局最低能量值为-2544.3431，全局最优结构为类笼型结构。由于异构体各原子之间满足一定的约束关系，因此预测得到的金团簇Au32最优结构较为稳定。
针对问题 3，传统模型无法很好地表征硼团簇的能量，为此我们采用机器学习方法对硼团簇能量及空间结构间的关系进行学习，提出了基于库伦矩阵本征谱的能量预测模型。结合 AWPSO 算法改进其适应度函数，设计了基于神经网络适应度函数的 AWPSO全局寻优算法。结果表明，硼团簇B45− 的最低能量值为-114305.0794，全局最优结构为带有一定曲率的类平面结构，中间有一六边形空位。
针对问题 4，我们改进了问题 2 的寻优方法，结合库伦矩阵本征谱，提出了多条件约束的 NN-AWPSO 模型。结果表明，硼团簇B40− 的最低能量值为-101954.52，全局最优结构为类平面结构，中间存在两个六边形空位，每个原子均可与临近原子构成类正六边形结构，相较于B45− 该结构对称性与稳定性较高。
最后，我们对模型进行了灵敏度分析不断改变学习因子的数值，发现模型预测出的能量值波动幅度较低，具有较强的稳定性，证明了建立模型的可靠性、有效性及鲁棒性，并进行了模型的评价与推广。
模型假设：
为便于建立团簇的数学模型与最优构型的求解，在建模过程中，进行如下假设：
(1) 所给附件中团簇势能未注明单位，为统一所有问题的分析与求解，将金团簇与硼团簇的势能值单位转换为附件所给对应势能值的相同量纲。
(2) 只考虑原子个数、相对位置对于团簇结构势能的影响，不考虑原子之间的角度对势能的影响。
(3) 假设对于同一种元素构成的团簇，其势能函数中与势能计算有关的常值参数保持不变。
(4) 在计算团簇势能时，假设团簇外部环境对该团簇的势能没有影响，即团簇处在一个孤立的空间内。
问题分析：
问题 1 的思路分析
本题中，要求根据所给附件中 1000 个金团簇𝐴𝑢20的结构及对应能量，建立金团簇的能量预测模型。对于金属团簇，物理化学领域的学者们提出了多种能够较好地描述团簇势能关系的模型，如𝐿𝑒𝑛𝑛𝑎𝑟𝑑 − 𝐽𝑜𝑛𝑒𝑠势能模型[1]、𝑀𝑜𝑟𝑠𝑒势能模型[2]、𝐺𝑢𝑝𝑡𝑎势能模型[3]、𝑆𝑢𝑡𝑡𝑜𝑛 − 𝐶ℎ𝑒𝑛势能模型[4]等，能将这一物理化学问题转化成纯粹的数学问题。因此，应选择合适的模型，并加入变异因子，采用机器学习非线性回归的方法，得到参数值，构成金团簇能量预测模型。针对预测金团簇𝐴𝑢20的全局最优结构的任务要求，考虑到粒子群优化算法(Particle Swarm optimization, PSO)作为一种具有全局寻优能力的群智能算法，本文对传统 PSO 算法中的惯性权重w进行改进，使其在探索阶段很大，而在开发阶段减小，形成自适应权重粒子群优化算法(Adaptive Weight Particle Swarm optimization, AWPSO)。使用 AWPSO算法预测金团簇𝐴𝑢20的全局最优结构，并对预测结果的几何结构和特点进行描述。
问题 2 的思路分析
本题要求在问题 1 的基础上，设计能够产生金团簇不同结构异构体的算法。考虑到团簇中各个原子的绝对位置坐标信息复杂，无法直观地表现原子间的几何关系，因此应建立相对距离矩阵(Relative Distance Matrix, RDM)，来描述团簇中各原子间的空间位置关系。不同的相对距离矩阵，即代表金团簇的不同结构异构体。另外，通过查阅文献可知，团簇的全局稳定结构也叫基态结构，一般情况下，结构的对称性越强，其势能也越低，结构也越稳定。并且，团簇中两两原子间的相对距离也有范围要求，相对距离过大，原子间的作用力微弱，无法形成稳定结构；相对距离过小，原子间的斥力变大，团簇势能增加，无法构成基态结构。根据对称性及相对距离的约束条件，本题中我们提出了基于相对距离矩阵的多条件约束异构生成方法，能够产生金团簇不同结构异构体，结合问题1 中的 AWPSO 算法，预测出金团簇Au32的全局最优结构，并描述其几何结构特点。
问题 3 的思路分析
本题中，研究对象从金属金团簇转变为非金属硼团簇，二者的区别在于，对于非金属团簇，现有的传统势能模型无法很好的表征非金属团簇的势能，因而需要利用神经网络等机器学习的方法，探索硼团簇分子势能和其空间结构间的关系。库伦矩阵(Coulomb Matrix)能把由团簇结构的原子坐标信息和原子核电荷数构成的数据转化为n *n维的数值矩阵，其中n为团簇中的原子个数。提取库伦矩阵中的本征谱，将得到的n维本征谱向量与团簇能量值之间进行机器学习神经网络训练，从而构建基于库伦矩阵本征谱的硼团簇能量预测模型。结合问题 1 中使用的 AWPSO 算法，将其中的适应度函数改进为训练好的神经网络，形成基于神经网络适应度函数的 AWPSO 寻优模型(Neutral Network Adaptive Function-AWPSO, NN-AWPSO)。使用 NN-AWPSO 模型预测硼团簇B45− 的全局最优结构，并对预测结果的几何结构和特点进行描述。
问题 4 的思路分析
结合问题 2，本题考虑其对称性限制，将不同的库伦矩阵本征谱向量作为硼团簇不同结构的信息，并约束本征值的最大值最小值，从而生成硼团簇的异构体。利用问题 3中提出的 NN-AWPSO 算法，预测出硼团簇B40− 的全局最优结构，并描述其几何结构特点。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：(代码和文档not free)
The actual procedure is shown in the screenshot
import numpy as np
import mpl_toolkits.mplot3d
import matplotlib.pyplot as plt
import pandas as pd
plt.figure(figsize=(24, 10))
df = pd.read_excel("0.xlsx")
x=df['x']
y=df['y']
z=df['z']
ax=plt.subplot(111,projection='3d')
for i in range(len(x)):
for j in range(len(y)):
ax.plot((x[i],x[j]),(y[i],y[j]),(z[i],z[j]),color='pink')
for i in range(len(x)):
ax.text(x[i],y[i],z[i],i,color='blue')
ax.set_zlabel('z')
ax.set_ylabel('y')
ax.set_xlabel('x')
plt.show()
#ML method Au20
#present path file name
import os
import pandas as pd
import re
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split
#find file name
file_name = os.listdir()
for f_name in file_name:
#judge file is not txt file
if f_name.endswith('.txt'):
# read data
with open(f_name,'r',encoding='utf-8') as f:
data = f.readlines()
# data slice get number type data
data = data[2:]
# get number use re to get
number_str_list = [re.findall(r"[- ]\d+\.?\d*",data[i].replace('\n','').strip(' ')) for i in range(20)]
#list store float data
number_float_list = []
#list store data label
number_target = []
# change str to float type
for item in number_str_list:
temp = []
count = 0
for n in item:
#change str data to float type
temp.append(float(n))
n = float(n)
#set data label
if n < 0:
count = count +1
#store data label
#number < 0 :number is 2 set label 0
if count > 1:
number_target.append(0)
#number < 0 :number is 1 set label 1
elif 0 < count < 2:
number_target.append(1)
##number < 0 :number is 0 set label 2
else:
number_target.append(2)
number_float_list.append(temp)
#split data
x_train,x_test, y_train, y_test =train_test_split(number_float_list,number_target,test_size=0.25, random_state=0)
#fit train data
svm = SVC().fit(x_train,y_train)
#use model to predict
predict = svm.predict(x_test)
#accuracy data set 1
acc = [1 for i in range(len(y_test)) if predict[i] == y_test[i]]
#cal accuracy
accuracy = len(acc) / len(y_test)
#output filename and accuracy
print(f_name,'--->',accuracy)
#Au20 ML
#import some tools
import os
import pandas as pd
import re
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split
#model train and predict
def modelSvm(number_float_list,number_target):
#split data
x_train,x_test, y_train, y_test =train_test_split(number_float_list,number_target,test_size=0.25, random_state=0)
#fit train data
svm = SVC().fit(x_train,y_train)
#use model to predict
predict = svm.predict(x_test)
#return predict label
return predict,y_test
#find file name
file_name_list = os.listdir()
for fname in file_name_list:
#judge file is not txt file
if fname.endswith('.txt'):
# read data
with open(fname,'r',encoding='utf-8') as f:
data = f.readlines()
# data slice get number type data
data = data[2:]
# get number use re to get
number_str_list = [re.findall(r"[- ]\d+\.?\d*",data[i].replace('\n','').strip(' ')) for i in range(45)]
#list store float data
data_float_list = []
#list store data label
data_target = []
# change str to float type
for item in number_str_list:
temp = []
count = 0
for n in item:
#change str data to float type
temp.append(float(n))
n = float(n)
#set data label
if n < 0:
count = count +1
#store data label
#number < 0 :number is 2 set label 2
if count > 1:
data_target.append(2)
#number < 0 :number is 1 set label 1
elif 0 < count < 2:
data_target.append(1)
##number < 0 :number is 0 set label 0
else:
data_target.append(0)
data_float_list.append(temp)
#use model to predict
predict,y_test = modelSvm(data_float_list,data_target)
#accuracy data set 1
acc = [1 for i in range(len(y_test)) if predict[i] == y_test[i]]
#cal accuracy
accuracy = len(acc) / len(y_test)
#output filename and accuracy
print(fname,'--->',accuracy)
752
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
