import random
import numpy as np
import sympy
arry11=np.zeros((20,20))
def func1_V2(data_distance): #计算V2
V1=np.zeros(20)
W=np.zeros(20)
a=np.zeros((20,20))
b=np.zeros((20,20))
c=np.zeros(20)
V2=np.zeros(20)
for i in range(20):
for j in range(20):
a[i][j]=0.2061*sympy.exp((-10.229) * (data_distance[i][j])-1)
V1[i]=V1[i]+a[i][j]
for i in range(20):
for j in range(20):
b[i][j]=1.79*sympy.exp((-2*4.036) * (data_distance[i][j])-1)
V2[i]=V2[i]+b[i][j]
return V2
def func_distance(arry11): #原子之间的距离公式
data_distance = np.zeros((20, 20))
for i in range(20):
for j in range(20 - i):
data_distance[i][j] = sympy.sqrt((arry11[i][0] - arry11[j][0]) ** 2 + (arry11[i][1] - arry11[j][1]) ** 2 + (arry11[i][2] - arry11[j][2]) ** 2)
return data_distance
def func_gradient(data_distance): #计算梯度
V2=func1_V2(data_distance)
r1 = np.zeros((20, 20))
r2 = np.zeros((20, 20))
r3 = np.zeros((20, 20))
for i in range(20):
for j in range(20):
r1[i][j] = 0.2061 * (-10.229) * sympy.exp((-10.229) * (data_distance[i][j] - 1))
r2[i][j] = 0.5 * V2[i] * 1.79 * (-2 * 4.036) * sympy.exp((-2 * 4.036) * (data_distance[i][j] - 1))
r3[i][j] = r1[i][j] - r2[i][j]
return r3
def func1_Gupta(data_distance): #计算Gupta势能
V1=np.zeros(20)
a=np.zeros((20,20))
b=np.zeros((20,20))
c=np.zeros(20)
V2=np.zeros(20)
V3 = np.zeros(20)
V_W=0
for i in range(20):
for j in range(20):
a[i][j]=0.2061*sympy.exp((-10.229) * (data_distance[i][j])-1)
V1[i]=V1[i]+a[i][j]
for i in range(20):
for j in range(20):
b[i][j]=1.79*sympy.exp((-2*4.036) * (data_distance[i][j])-1)
V2[i]=V2[i]+b[i][j]
for i in range(20):
V3[i]=sympy.sqrt(V2[i])
for i in range(20):
c[i]=V1[i]-V3[i]
for i in range(20):
V_W=V_W+c[i]
return V_W
def DFP(arry11): #局部极小化
a=sympy.symbols("a")
r=np.zeros((20,20))
r1=np.zeros((20, 20))
H=np.eye(20)
Nolmprove=0
data_distance1=np.zeros((20,20))
data_distance1=func_distance(arry11)
while True:
r1=func_gradient(data_distance1)
p=-np.dot(H,r1)
arry11 = arry11
for i in range(20):
for j in range(3):
arry11[i][j]=arry11[i][j]+a
alpha=func1_Gupta(arry11)
difyL_a = sympy.diff(alpha, a)
aa = sympy.solve([difyL_a], [a])
a = aa.get(a)
arry11_1=arry11
arry11=arry11+np.dot(a,p)
arry11_2=arry11
data_distance2=func_distance(arry11)
r2=func_gradient(data_distance2)
for i in range(20):
for j in range(20):
r[i][j]=r2[i][j]-r1[i][j]
arry11_3=arry11_2-arry11_1
H1=np.dot(arry11_3,np.transpose(arry11_3))/np.dot(np.transpose(arry11_3),r)-np.dot(H,np.dot(r,np.transpose(np.dot(H,r))))/np.dot(transpose(r),np.dot(H,r))
H=H+H1
Nolmprove=Nolmprove+1
if Nolmprove ==40:
return arry11
break
else:
continue
if __name__ == '__main__':
print("开始运行！")
arry11=np.zeros((20,20))
for i in range(20):
for j in range(20):
arry11[i][j]=random.uniform(-2,2)
data_distance_1 = np.zeros(20)
arry11=DFP(arry11)
arry11_1 = arry11
data_distance = func_distance(arry11)
Gupta1 = func1_Gupta(data_distance)
max_1 = 0
max_distance = 0
x = 0
y = 0
z = 0
x_list =np.zeros(20)
y_list =np.zeros(20)
z_list =np.zeros(20)
while True:
for i in range(20):
for j in range(3):
arry11[i][j] = arry11[i][j] + random.uniform(-0.5, 0.5);
for i in range(20):
x_list[i] = arry11[i][0]
x = x + x_list[i]
y_list[i] = arry11[i][1]
y = y + y_list[i]
z_list[i] = arry11[i][2]
z = z + z_list[i]
x = x / 20
y = y / 20
z = z / 20
xyz = np.zeros(3)
for i in range(20):
data_distance[i] = sympy.sqrt((arry11[i][0] - x) ** 2 + (arry11[i][1] - y) ** 2 + (arry11[i][2] - z) ** 2);
if data_distance >= max_distance:
max_distance = data_distance
xyz[0] = arry11[i][0]
xyz[1] = arry11[i][1]
xyz[2] = arry11[i][2]
a = i
arry11[a][0] = xyz[0] + random.uniform(-0.3, 0.3)
arry11[a][1] = xyz[1] + random.uniform(-0.3, 0.3)
arry11[a][2] = xyz[2] + random.uniform(-0.3, 0.3)
arry11 = DFP(arry11)
data_distance = func_distance(arry11)
Gupta2 = func1_Gupta(data_distance)
if Gupta2 <= Gupta1:
Gupta1 = Gupta2
arry11 = arry11
arry11_1 = arry11
max_1 = max_1 + 1
continue
else:
arry11 = arry11_1 # 回到初始
if max_1 == 30:
print(arry11)
break
else:
continue
DFP算法解法出现问题，运行出错，调用子程序时当返回值为数组时主程序无法读入，后续操作无法进行
求大佬给出有效代码
把82行到85行换成
arraySympy = np.array([])
for i in range(20):
for j in range(3):
arraySympy = np.append(arraySympy,arry11[i][j]+a)
print(arraySympy)