python 拟合 Lennard-Jones 对函数
python fit the Lennard-Jones pair function
我有一个双星系统(AB),并计算了10个结构的总能量。然后,我想使用 python 来使用 Lennard-Jones (LJ) 对函数拟合能量并获得它们相应的参数(epsilon 和 sigma)。 LJ函数是:
U(r)=4*epsilon [(sigma/r)^12 + (sigma/r)^6],
其中 r 是原子间距离。对于AB双星系统,它包括三种类型的相互作用:A-A、A-B和B-B,因此有六个拟合参数epsilon(A-A)、sigma(A-A)、epsilon(A-B)、sigma(A-B)、epsilon(B-B)和sigma( B-B)是必需的。换言之,总的U函数中有[U(A-A)、U(A-B)和U(B-B)]三部分。输入 r 变量是这样的列表:
r=[[[A,A],2.0], [[A,A],3.0],...,[[A,B],2.0], [[A,B],3.0],...,[[B,B],2.0],[[B,B],3.0] ...]
因此,我需要将 A-A、A-B 和 B-B 原子间距离分组到相应的函数 (U(A-A)、U(A-B) 和 U(B-B))。但是,我不知道如何正确编写函数和调用行:
def func(r,energy):
for i in range(len(r)):
if r[i][0]==[A,A]:
U0=U(A-A)
elif: r[i][0]=[A,B]:
U0=U(A-B)
else: U0=U(B-B)
U=U+U0
return U
popt, pcov = curve_fit(func, r, energy)
另一方面,我也想同时拟合力(势的一阶导数)?如何同时最小化能量和力拟合误差,最小化(能量,力),以获得能量和力的最佳拟合参数(epsilon,sigma)?非常感谢你。
我相信你想要的是一个 "atomic pairs" 的列表,每个都有 3 个值:(type1, type2, distance**6),其中 type1 和 type2 每个都有两个值之一 'a' 或 'b'。您 可以 从 "atomic coordinates" 的列表开始,(类型,x,y,z),但您可以随后将该列表缩减为 "atomic pairs" 的列表。
然后,您需要 6 个变量:epsAA、epsAB、epsBB、sigAA、sigAB、sigBB。
你的函数应该将一个数组作为输入,其中包含这 6 个变量和对数据的值,以及 return 一个能量数组(即每对的能量),可能类似于:
import numpy as np
def objective(params, pairsdata):
epsAA, epsAB, epsBB, sigAA, sigAB, sigBB = params
u = np.zeros(len(pairsdata))
for i, pairs in enumerate(pairsdata):
t1, t2, r6 = pairsdata
eps, sig6 = epsAB, sigAB # default to different types
if t1 == t2:
eps, sig6 = epsAA, sigAA
if t1 == 'b':
eps, sig6 = epsBB, sigBB
arg = sig6/r6
u[i] = 4 * eps * arg * (arg + 1)
return u
然后,您可以将其与 scipy.optimize.leastsq()
一起使用。
我想你可能想存储 r**6
(和 sigma**6
),然后平方,而不是每次都使用六次方和十二次方(这些 Lennard 和 Jones 的人在想什么? 计算机具有无限的精度?;))。
值得考虑如何避免 for 循环,但你并没有说你希望它快,只是正确 ;)。
如果你想包括一些力的度量,我认为你需要通过有限差分法计算,然后通过 concatenating 能量数组和力数组(即 return 数组的值是对的两倍)。
感谢纽维尔的回答。对于更简单的情况,我编写了一个简短的 python 程序。它可能会澄清我的问题。
import numpy as np
import scipy.optimize
def func(p,x,ydata):
f=4*p[0]*( (p[1]/x)**12 - (p[1]/x)**6 )
return f
dis=[[3.45454545455,3.63636363636,4.0],[3.54545454545,4.18181818182,5.0],[3.81818181818,4.45454545455,4.90909090909], [3.72727272727,4.36363636364,5.36363636364],[3.90909090909,5.27272727273,6.0]]
ene=[-1.3, -1.4, -2.0, -2.2, -1.3]
xdata=np.array(dis)
ydata=np.array(ene)
p0=[0.5,3.0]
print scipy.optimize.leastsq(func, p0, args=(xdata,ydata))
我有一个双星系统(AB),并计算了10个结构的总能量。然后,我想使用 python 来使用 Lennard-Jones (LJ) 对函数拟合能量并获得它们相应的参数(epsilon 和 sigma)。 LJ函数是:
U(r)=4*epsilon [(sigma/r)^12 + (sigma/r)^6],
其中 r 是原子间距离。对于AB双星系统,它包括三种类型的相互作用:A-A、A-B和B-B,因此有六个拟合参数epsilon(A-A)、sigma(A-A)、epsilon(A-B)、sigma(A-B)、epsilon(B-B)和sigma( B-B)是必需的。换言之,总的U函数中有[U(A-A)、U(A-B)和U(B-B)]三部分。输入 r 变量是这样的列表:
r=[[[A,A],2.0], [[A,A],3.0],...,[[A,B],2.0], [[A,B],3.0],...,[[B,B],2.0],[[B,B],3.0] ...]
因此,我需要将 A-A、A-B 和 B-B 原子间距离分组到相应的函数 (U(A-A)、U(A-B) 和 U(B-B))。但是,我不知道如何正确编写函数和调用行:
def func(r,energy):
for i in range(len(r)):
if r[i][0]==[A,A]:
U0=U(A-A)
elif: r[i][0]=[A,B]:
U0=U(A-B)
else: U0=U(B-B)
U=U+U0
return U
popt, pcov = curve_fit(func, r, energy)
另一方面,我也想同时拟合力(势的一阶导数)?如何同时最小化能量和力拟合误差,最小化(能量,力),以获得能量和力的最佳拟合参数(epsilon,sigma)?非常感谢你。
我相信你想要的是一个 "atomic pairs" 的列表,每个都有 3 个值:(type1, type2, distance**6),其中 type1 和 type2 每个都有两个值之一 'a' 或 'b'。您 可以 从 "atomic coordinates" 的列表开始,(类型,x,y,z),但您可以随后将该列表缩减为 "atomic pairs" 的列表。
然后,您需要 6 个变量:epsAA、epsAB、epsBB、sigAA、sigAB、sigBB。
你的函数应该将一个数组作为输入,其中包含这 6 个变量和对数据的值,以及 return 一个能量数组(即每对的能量),可能类似于:
import numpy as np
def objective(params, pairsdata):
epsAA, epsAB, epsBB, sigAA, sigAB, sigBB = params
u = np.zeros(len(pairsdata))
for i, pairs in enumerate(pairsdata):
t1, t2, r6 = pairsdata
eps, sig6 = epsAB, sigAB # default to different types
if t1 == t2:
eps, sig6 = epsAA, sigAA
if t1 == 'b':
eps, sig6 = epsBB, sigBB
arg = sig6/r6
u[i] = 4 * eps * arg * (arg + 1)
return u
然后,您可以将其与 scipy.optimize.leastsq()
一起使用。
我想你可能想存储 r**6
(和 sigma**6
),然后平方,而不是每次都使用六次方和十二次方(这些 Lennard 和 Jones 的人在想什么? 计算机具有无限的精度?;))。
值得考虑如何避免 for 循环,但你并没有说你希望它快,只是正确 ;)。
如果你想包括一些力的度量,我认为你需要通过有限差分法计算,然后通过 concatenating 能量数组和力数组(即 return 数组的值是对的两倍)。
感谢纽维尔的回答。对于更简单的情况,我编写了一个简短的 python 程序。它可能会澄清我的问题。
import numpy as np
import scipy.optimize
def func(p,x,ydata):
f=4*p[0]*( (p[1]/x)**12 - (p[1]/x)**6 )
return f
dis=[[3.45454545455,3.63636363636,4.0],[3.54545454545,4.18181818182,5.0],[3.81818181818,4.45454545455,4.90909090909], [3.72727272727,4.36363636364,5.36363636364],[3.90909090909,5.27272727273,6.0]]
ene=[-1.3, -1.4, -2.0, -2.2, -1.3]
xdata=np.array(dis)
ydata=np.array(ene)
p0=[0.5,3.0]
print scipy.optimize.leastsq(func, p0, args=(xdata,ydata))