Вычисление параметров для молекулярной механики

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize 
import subprocess
In [3]:
!orca

This program requires the name of a parameterfile as argument
For example ORCA TEST.INP

orca работает

Длина связи

Изменяем длину связи от 1.2 до 1.7 ангстрем с шагом 0.02.

In [4]:
def run_orca():
    with subprocess.Popen("orca orca.inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
        out=p.communicate()[0].decode()
    print(out)
    for line in out.splitlines():
        if 'FINAL SINGLE POINT' in line:
            return float(line.strip().split()[4])
    return None
In [ ]:
bond_E = []
for length in np.round(np.arange(1.2,1.7,0.02),2):
    inp = f'''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 {length} 0 0 
    H 1 2 0 1.08439 111.200 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 180
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    '''
    with open('orca.inp','w') as newfile:
        newfile.write(inp)
    energy = run_orca()
    bond_E.append(energy)

Посмотрим, как выглядит зависимость.

In [6]:
f, ax = plt.subplots(1,1, figsize=(7, 7))
plt.scatter(np.round(np.arange(1.2,1.7,0.02),2), bond_E)
Out[6]:
<matplotlib.collections.PathCollection at 0x7f357dd24b00>
In [15]:
#function is  f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,-79, 1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc,
                               p0[:],
                               args=(np.round(np.arange(1.2,1.7,0.02),2), bond_E))
print("Optimized params:", p1)

#Plot it
f, ax = plt.subplots(1,1, figsize=(7, 7))
plt.plot(np.round(np.arange(1.2,1.7,0.02),2),
         bond_E,
         "ro",
         np.round(np.arange(1.2,1.7,0.02),2),
         fitfunc(p1, np.round(np.arange(1.2,1.7,0.02),2)))

plt.xlim(1.1,1.8)
plt.show()

Optimized params: [  0.9917356    1.53476476 -79.20163191]

Зависимость не очень хорошо аппроксимируется функцией. Попробуем выкинуть несколько точек.

In [22]:
#function is  f(x)=k(b-x)^2 + a
x = np.round(np.arange(1.2,1.7,0.02),2)[9:]
y = bond_E[9:]
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,-79, 1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc,
                               p0[:],args=(x,y))
print("Optimized params:", p1)

#Plot it
f, ax = plt.subplots(1,1, figsize=(7, 7))
plt.plot(x,
         y,
         "ro",
         x,
         fitfunc(p1, x))

plt.xlim(1.1,1.8)
plt.show()

Optimized params: [  0.59809713   1.54252522 -79.19772785]

Выкинули первые 9 точек, и функция стала фиттиться гораздо лучше.

Валентный угол

Изменяем угол HCC от 109.2 до 113.4 с шагом 0.1.

In [ ]:
angle_E = []
for angle in np.round(np.arange(109.2,113.4,0.1),1):
    inp = f'''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 1.52986 0 0 
    H 1 2 0 1.08439 {angle} 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 180
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    '''
    with open('orca.inp','w') as newfile:
        newfile.write(inp)
    energy = run_orca()
    angle_E.append(energy)
In [24]:
f, ax = plt.subplots(1,1, figsize=(7, 7))
plt.scatter(np.round(np.arange(109.2,113.4,0.1),1), angle_E)
Out[24]:
<matplotlib.collections.PathCollection at 0x7f357c202c18>
In [25]:
#function is  f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,-79, 1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc,
                               p0[:],
                               args=(np.round(np.arange(109.2,113.4,0.1),1), angle_E))
print("Optimized params:", p1)

#Plot it
f, ax = plt.subplots(1,1, figsize=(7, 7))
plt.plot(np.round(np.arange(109.2,113.4,0.1),1),
         angle_E,
         "ro",
         np.round(np.arange(109.2,113.4,0.1),1),
         fitfunc(p1, np.round(np.arange(109.2,113.4,0.1),1)))

plt.xlim(108,114)
plt.show()

Optimized params: [ 4.05989499e-05  1.11195346e+02 -7.91975724e+01]

Хорошо зафиттилось. Точки выкидывать не будем.

Торсионный угол

Изменяем торсионный угол от -180 до 180 с шагом 12.

In [ ]:
torsion_E = []
for torsion in np.round(np.arange(-180,181,12),1):
    inp = f'''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 1.52986 0 0 
    H 1 2 0 1.08439 111.200 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 {torsion}
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    '''
    with open('orca.inp','w') as newfile:
        newfile.write(inp)
    energy = run_orca()
    torsion_E.append(energy)
In [27]:
#function is  f(x)=k/2 * (1+cos(n*x - g)) + a

fitfunc = lambda p, x: p[0] / 2 * (1 + np.cos(np.radians(p[1] * np.pi * x - p[2]))) + p[3]# Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,1, 1,1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc,
                               p0[:],
                               args=(np.round(np.arange(-180,
                                                        181,
                                                        12),1), torsion_E))
print("Optimized params:", p1)

#Plot it
f, ax = plt.subplots(1,1, figsize=(7, 7))
plt.plot(np.round(np.arange(-180,181,12),1),
         torsion_E,
         "ro",
         np.round(np.arange(-180,181,12),1),
         fitfunc(p1, np.round(np.arange(-180,181,12),1)))

plt.show()
    

Optimized params: [ 4.58410226e-03  9.54986108e-01  7.94011482e-05 -7.91975763e+01]

Зафиттилось тоже хорошо.

Длина связи

Изменяем длину связи от 0.1 до 5 с шагом 0.1.

In [ ]:
bond_E2 = []
for length in np.round(np.arange(0.1,5,0.1),2):
    inp = f'''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 {length} 0 0 
    H 1 2 0 1.08439 111.200 0
    H 1 2 3 1.08439 111.200 120
    H 1 2 3 1.08439 111.200 -120
    H 2 1 3 1.08439 111.200 180
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    '''
    with open('orca.inp','w') as newfile:
        newfile.write(inp)
    energy = run_orca()
    bond_E2.append(energy)
In [32]:
f, ax = plt.subplots(1,1, figsize=(7, 7))
plt.scatter(np.round(np.arange(0.1,5,0.1),2), bond_E2)
Out[32]:
<matplotlib.collections.PathCollection at 0x7f35767ee048>

Нужно выкинуть начальные точки, потому что из-за их очень низкого значения (непонятно, почему такого низкого) не видно, как выглядит функция. А функция должна быть похожа на потенциал Морзе, потому что именно им в данном случае ее нужно аппроксимировать.

In [35]:
# morze potential

fitfunc = lambda p, x: p[0] * (1 - np.exp(- p[1] * (x - p[2]))) ** 2 + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
x = np.round(np.arange(0.1,5,0.1),2)[4:]
y = bond_E2[4:]

p0 = [1,1,1,1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x,y))

print("Optimized params:", p1)

#Plot it
f, ax = plt.subplots(1,1, figsize=(7, 7))
plt.plot(x, y, "ro", x, fitfunc(p1, x))
plt.show()
    

Optimized params: [  0.30555985   1.90108474   1.40990572 -79.21751735]

Гораздо лучше без первых четырёх точек.

In [1]:
from IPython.core.display import display, HTML

html_custom = '<h4>Файлы, полученные выше, лежат в этой <a href="https://kodomo.fbb.msu.ru/~agayeva_zarifa/term8/pr6_files">папке</a>.</h4>'

display(HTML(html_custom))

Файлы, полученные выше, лежат в этой папке.