S.M. Ikhdair, et al., J. Math. Chem., 41 (2007), p. 329.
https://link.springer.com/content/pdf/10.1007%2Fs10910-006-9080-2.pdf
In [1]:
import symengine as se
one = lambda: se.Rational(1, 1)
nL = one()*0
nA = one()*1
ndelta = one()*10/100
nhbar, nm = one()*1, one()*1
In [2]:
En0 = lambda n: -nm*nA**2/(2*nhbar**2*(n+nL+1)**2)
E01 = -nhbar**4*(nL+1)**2*(nL+2)*(2*nL+3)*ndelta**3/(6*nA*nm**2)
E021 = nhbar**6*(nL+1)**3*(nL+2)*(2*nL+3)*(2*nL+5)*ndelta**4/(24*nA**2*nm**3)
E022 = -nhbar**10*(nL+1)**6*(nL+2)*(2*nL+3)*(8*nL**2+37*nL+43)*ndelta**6/(72*nA**4*nm**5)
E02 = E021 + E022
print((En0(0) + nA*ndelta + E01 + E02)*1.0)
In [ ]: