0

说到这，我很绿色sympy，我不知道如何以格式良好的方式生成输出。现在我已经计算了我的潜在函数的 Hessian 矩阵：

V = 1/2*kOH*(r1)**2 +1/2*kOH*(r2)**2 +1/2*kHH*(r3)**2

三个谐振子项的一般形式为：

1/2*k*r**2.

所有变量都是正数和实数。

对我来说问题是，当我打印我的矩阵时，条目还没有解决，只以功能方式显示。我希望在已经执行偏导数之后将条目放在表格中，而不仅仅是显示在矩阵中的每个点需要执行哪些推导。

def Hessian():

    '''
    sympy calc of hessian Matrix H for IR normal modes analysis
    from a potential V.

    Must be multiplicable with 9x9 matrix (somehow) in 
    the equation: F = M**(-1/2) * H * M**(-1/2)
    Here, F is the mass weighted Hessian, whose Eigenvalues
    contain the frequencies of the normal modes of water.
    M comes from the multiplication of the 3N-Dimensional
    mass-vector m with a 3N-dimensional identity matrix:
    M = m*I, I.shape = 3*N, 3*N, N = number of atoms in water. 
    '''

    kOH, kHH, r1, r2, r3 = sy.symbols('kOH kHH r1 r2 r3', real=True, positive=True)

    V = sy.Function('V')(1/2*kOH*(r1)**2 +1/2*kOH*(r2)**2 +1/2*kHH*(r3)**2)

    f = sy.hessian(V,[r1, r2, r3])

    sy.pprint(f)


Hessian()

附加：这实际上不是事物计算方面的一部分，因此也不是问题的一部分，但是如果有人在事物的科学方面知道他们的东西：你能告诉我如何 a (3,3)潜在依赖于三个距离的 Hessian 矩阵应该乘以 (9,9) 质量矩阵？如果您有兴趣，函数的注释包含科学背景。

4

1 回答 1

1

您遇到的基本问题是：

In [39]: f = Function('f')                                                                                                        

In [40]: f(x)                                                                                                                     
Out[40]: f(x)

In [41]: f(x).diff(x)                                                                                                             
Out[41]: 
d       
──(f(x))
dx      

In [42]: f(x).diff(x).subs(x, 2*y)                                                                                                
Out[42]: 
⎛d       ⎞│     
⎜──(f(x))⎟│     
⎝dx      ⎠│x=2⋅y

理想情况下，SymPy 会将最后一个结果表示为类似的东西，f'(2y)但 SymPy 没有办法直接表示这样的对象。理想情况下，会有一个微分运算符D，这样 sayD(f)(x)将与 相同f(x).diff(x)。有了它，您可以将D(f)(2*y)其表示为当然可以显示为f'(2y).

当然，如果你用一个函数代替f这里，那么导数可以评估：

In [45]: f(x).diff(x).subs(x, 2*y).subs(f, Lambda(t, t**3))                                                                       
Out[45]: 
⎛d ⎛ 3⎞⎞│     
⎜──⎝x ⎠⎟│     
⎝dx    ⎠│x=2⋅y

In [46]: _.doit()                                                                                                                 
Out[46]: 
    2
12⋅y

要回答您的另一个问题，显然您不能将 9x9 矩阵和 3x3 矩阵相乘。您的等式F意味着H和M都是正方形且大小相同。您的质量矩阵实际上只是 3x3，或者您的潜在函数实际上是 9 个坐标的函数。假设r1是原子 1 和原子 2 之间的距离，那么r1 = sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)在这种情况下，你应该计算你的 Hessian wrtx1等而不是r1.

于 2019-09-18T20:20:48.950 回答