Python Script For Solving Mp2 Equations with Solution is the subject that will revolve around this post. This article solves linear and nonlinear equations with Python programming. The result of linear equations is through the matrix operations while sets of nonlinear equations make it necessary for a solver numerically find a solution.
2x1+3x2+4x3=5
5x1+6x2+7x3=8
8x1+10x2=12
First of all import NumPy and assign value as a set.
import numpy as np
A = np.array([[2, 3, 4],
[5, 6, 7],
[8, 10, 0]])
y = np.array([5, 8, -12])
x = np.linalg.solve(A, y)
print(x)
The output of the 2x1+3x2+4x3=5
5x1+6x2+7x3=8
8x1+10x2=12 Is following .
[ 0.16666667 -1.33333333 2.16666667]
Find the 2x1+3x2+4x3=5
5x1+6x2+7x3=8
8x1+10x2=12 Using the matrix inversion approach
A_inv = np.linalg.inv(A)
x = np.dot(A_inv, y)
print(x)The out-of-code is following
[ 0.16666667 -1.33333333 2.16666667]
Find the LL and UU for the above matrix A.
from scipy.linalg import lu
P, L, U = lu(A)
print(‘P:\n’, P)
print(‘L:\n’, L)
print(‘U:\n’, U)
print(‘LU:\n’,np.dot(L, U))
P:
[[0. 1. 0.]
[0. 0. 1.]
[1. 0. 0.]]
L:
[[ 1. 0. 0. ]
[ 0.25 1. 0. ]
[ 0.625 -0.5 1. ]]
U:
[[ 8. 10. 0. ]
[ 0. 0.5 4. ]
[ 0. 0. 9. ]]
LU:
[[ 8. 10. 0.]
[ 2. 3. 4.]
[ 5. 6. 7.]]
Find PP and AA and see what’s the effect of the permutation matrix on AA.
The LL and UU we get are different from the ones and got in the last section by hand. And also see there is a permutation matrix PP that occurs after a period of absence by the lu function. This permutation matrix record how we change the order of the equations for easier calculation purposes (if the first element in the first row is zero, it cannot be the pivot equation, since you cannot turn the first elements in other rows to zero. We need to switch the order of the equations to get an add new pivot equation).
Also Read: Division Operator in Python – A Complete Guide
If you multiply PP with AA, you will see that this is each of several possible ways in which a set or number of things can be ordered matrix reverse the order of the equations for this case.
print(np.dot(P, A))
[[ 5. 6. 7.]
[ 8. 10. 0.]
[ 2. 3. 4.]]
import numpy as np
A = np.array([ [3,-9], [2,4] ])
b = np.array([-42,2])
z = np.linalg.solve(A,b)
print(z)
M = np.array([ [1,-2,-1], [2,2,-1], [-1,-1,2] ])
c = np.array([6,1,1])
y = np.linalg.solve(M,c)
print(y)
[-5. 3.]
[ 3. -2. 1.]
import numpy as np
from scipy.optimize import fsolve
def myFunction(z):
x = z[0]
y = z[1]
w = z[2]
F = np.empty((3))
F[0] = x**2+y**2-20
F[1] = y – x**2
F[2] = w + 5 – x*y
return F
zGuess = np.array([1,1,1])
z = fsolve(myFunction,zGuess)
print(z)
from gekko import GEKKO
m = GEKKO()
x,y,w = [m.Var(1) for i in range(3)]
m.Equations([x**2+y**2==20,\
y-x**2==0,\
w+5-x*y==0])
m.solve(disp=False)
print(x.value,y.value,w.value)
import sympy as sym
sym.init_printing()
x,y,z = sym.symbols(‘x,y,z’)
c1 = sym.Symbol(‘c1’)
f = sym.Eq(2*x**2+y+z,1)
g = sym.Eq(x+2*y+z,c1)
h = sym.Eq(-2*x+y,-z)
sym.solve([f,g,h],(x,y,z))
Python Script For Solving Mp2 Equations with Complete Solution
from __future__ import division
import math
import numpy as np
def eint(a,b,c,d):
if a > b: ab = a*(a+1)/2 + b
else: ab = b*(b+1)/2 + a
if c > d: cd = c*(c+1)/2 + d
else: cd = d*(d+1)/2 + c
if ab > cd: abcd = ab*(ab+1)/2 + cd
else: abcd = cd*(cd+1)/2 + ab
return abcd
def teimo(a,b,c,d):
return ttmo.get(eint(a,b,c,d),0.0e0)
Nelec = 2 # we have 2 electrons in HeH+
dim = 2 # we have two spatial basis functions in STO-3G
E = [-1.52378656, -0.26763148]
ttmo = {5.0: 0.94542695583037617, 12.0: 0.17535895381500544, 14.0: 0.12682234020148653, 17.0: 0.59855327701641903, 19.0: -0.056821143621433257, 20.0: 0.74715464784363106}
dim = dim*2
ints=np.zeros((dim,dim,dim,dim))
for p in range(1,dim+1):
for q in range(1,dim+1):
for r in range(1,dim+1):
for s in range(1,dim+1):
value1 = teimo((p+1)//2,(r+1)//2,(q+1)//2,(s+1)//2) * (p%2 == r%2) * (q%2 == s%2)
value2 = teimo((p+1)//2,(s+1)//2,(q+1)//2,(r+1)//2) * (p%2 == s%2) * (q%2 == r%2)
ints[p-1,q-1,r-1,s-1] = value1 – value2
fs = np.zeros((dim))
for i in range(0,dim):
fs[i] = E[i//2]
#fs = np.diag(fs)
EMP2 = 0.0
for i in range(0,Nelec):
for j in range(0,Nelec):
for a in range(Nelec,dim):
for b in range(Nelec,dim):
EMP2 += 0.25*ints[i,j,a,b]*ints[i,j,a,b]/(fs[i] +fs[j] -fs[a] – fs[b])
print(“E(MP2) Correlation Energy = “, EMP2, ” Hartrees”)
E(MP2) Correlation Energy = -0.0064020383431405556
Conclusion:
The result of linear equations is through the matrix operations while sets of nonlinear equations make necessary a solver numerically find a solution. The MP2 module supports spin-restricted orbitals and the CholeskyLinalgFactory and DenseLinalgFactory.