Python适合Lennard-Jones对函数

Question

我有一个双星系统，计算了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变量如下所示：

^{pr2}$

因此，我需要将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）？非常感谢你。在

Answer 1

我相信你想要的是一个“原子对”列表，每个原子对有3个值：（type1，type2，distance**6），其中type1和type2都有两个值“a”或“b”中的一个。你可以从一个“原子坐标”（type，x，y，z）列表开始，但是你可以把这个列表减少到“原子对”的列表中。在

然后，你需要6个变量：epsAA，epsAB，epsBB，sigAA，sigAB，sigBB。在

您的函数应该将包含这6个变量值的数组和成对数据作为输入，并返回一个能量数组（即，每对的能量），例如：

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），并将其平方，而不是每次都使用第六次和第十二次幂（这些伦纳德和琼斯的人在想什么？计算机有无限的精确性？;) ). 在

值得考虑如何避免for循环，但您并没有说您希望它更快，只是正确的；）。在

如果你想包括一些力的测量，我想你需要用一个有限差分法来计算，然后用U和dU/dr来连接能量数组和力的数组（即返回一个值是对的两倍的数组）。在

Answer 2

感谢纽维尔的回答。为了简单起见，我编写了一个简短的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))