I'm very green when it comes to sympy and I don't know how to produce output in a well formatted fashion. Right now I've computed the Hessian matrix of my potential function:
V = 1/2*kOH*(r1)**2 +1/2*kOH*(r2)**2 +1/2*kHH*(r3)**2
of three harmonic oscillator terms with the gereral form of:
1/2*k*r**2.
All variables are positive and real.
The problem for me is, that when I print my matrix, the entries are not yet solved and only displayed in a functional fashion. I'd like the entries to be in the form after the partial derivatives have already been preformed, not just displaying what derivations need to be preformed at each point in the matrix.
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()
Additional: This isn't really part of the computing side of things and so not really part of the question, but if anyone knows their stuff when it comes to the science side of things: could you tell me how a (3,3) Hessian of the potential dependend on three distances is supposed to be multiplied with a (9,9) Mass matrix? The commentary of the function contains the scientific background, if you're interested.
The basic problem you have is this:
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
Ideally SymPy would express that last result as something like f'(2y) but there SymPy doesn't have a way to represent such an object directly. Ideally there would be a differential operator D so that say D(f)(x) would be the same as f(x).diff(x). With that you could represent this as D(f)(2*y) which could of course display as f'(2y).
Of course if you substitute a function for f here then the derivative can evaluate:
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
To answer your other question, clearly you can't multiply a 9x9 matrix and a 3x3 matrix. Your equation for F implies that both H and M are square and of the same size. Either your mass matrix is really just 3x3 or your potential function is really a function of 9 coordinates. Assuming r1 is the distance between say atom 1 and atom 2 then maybe r1 = sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2) in which case you should compute your Hessian wrt x1 etc rather than r1.