python 拟合 Lennard-Jones 对函数

Question

我有一个双星系统（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）？非常感谢你。

Answer 1

我相信你想要的是一个 "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 数组的值是对的两倍）。

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))

python 拟合 Lennard-Jones 对函数

python fit the Lennard-Jones pair function

python

function

curve-fitting