Наша цель состоит в том, что бы рассчитать 20 разных длинн связи с шагом 0.02 ангстрема для расчёта энергии в ORCA с разными значениями по длине одной из связей (C-C)

In [7]:
%matplotlib inline
import subprocess
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize
In [4]:
#Для расчета разных значений длин связи испортируем функцию orca
def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    with subprocess.Popen("/srv/databases/orca/orca orca.inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
        out=p.communicate()[0].decode()
    for l in out.splitlines():
        if 'FINAL SINGLE POINT' in l:
            return float(l.strip().split()[4])
    return None
In [2]:
inp = '''!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 180
H 2 1 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
In [5]:
run_orca(inp)
Out[5]:
-79.081531418112

Сделаем шаблоны

In [3]:
bond = '''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0 
C 1 0 0 {:.5f} 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 5 1.08439 111.200 120
H 2 1 5 1.08439 111.200 -120
*
'''
angle = '''!MP2 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 {:.3f} 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
*'''
tors = '''!MP2 6-311G**
* 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 {:d}
H 2 1 6 1.08439 111.200 120
H 2 1 6 1.08439 111.200 -120
*'''
In [8]:
#сделаем массив с шагом по длине связи
bonds = np.arange(1.3, 1.7, 0.02)
bonds
Out[8]:
array([1.3 , 1.32, 1.34, 1.36, 1.38, 1.4 , 1.42, 1.44, 1.46, 1.48, 1.5 ,
       1.52, 1.54, 1.56, 1.58, 1.6 , 1.62, 1.64, 1.66, 1.68])
In [9]:
bond_energies = np.array([run_orca(bond.format(value))
                                 for value in bonds])
In [10]:
bond_energies
Out[10]:
array([-79.0321825 , -79.04174446, -79.04998776, -79.05703088,
       -79.06298212, -79.06794047, -79.07199644, -79.07523272,
       -79.0777249 , -79.07954205, -79.08074728, -79.08139822,
       -79.08154754, -79.0812433 , -79.08052939, -79.07944588,
       -79.07802931, -79.07631301, -79.07432738, -79.07210013])

Отрисуем зависимость длины связи от энергии

In [23]:
#функция  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, 1, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(bonds, bond_energies))
print("Optimized params:", p1)

#Plot it
plt.plot(bonds, bond_energies, 'ro',color='blue', alpha=0.5)
plt.plot(bonds, fitfunc(p1, bonds), "r-", c='red')
#plt.xlim(1,2)
plt.show()

Optimized params: [  0.77459065   1.54783147 -79.08299494]

Orca рассчитал связь в 1.54783147 ангстрем. В статье о разработке поля GAFF указана длина C-C 1.526. Разница довольно заметная, но тут вклад еще вносит неточность в фите - по результатам минимум приходится на меньшее значение, которое ближе к статье.

In [25]:
#сделаем массив с шагом по величине валентного угла HCС
angles = np.arange(109.2,113.3,0.2)
angles
Out[25]:
array([109.2, 109.4, 109.6, 109.8, 110. , 110.2, 110.4, 110.6, 110.8,
       111. , 111.2, 111.4, 111.6, 111.8, 112. , 112.2, 112.4, 112.6,
       112.8, 113. , 113.2])
In [26]:
angle_energies = np.array([run_orca(angle.format(value))
                                 for value in angles])
In [27]:
angle_energies
Out[27]:
array([-79.49175077, -79.49177525, -79.49179754, -79.49181698,
       -79.4918335 , -79.49184711, -79.49185781, -79.4918656 ,
       -79.49187049, -79.49187247, -79.49187155, -79.49186772,
       -79.49186101, -79.49185139, -79.49183889, -79.4918235 ,
       -79.49180521, -79.49178405, -79.49176   , -79.49173307,
       -79.49170326])

Отрисуем зависимость величины угла от энергии

In [29]:
#функция  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, 1, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(angles, angle_energies))
print("Optimized params:", p1)

#Plot it
plt.plot(angles, angle_energies, 'ro',color='blue', alpha=0.5)
plt.plot(angles, fitfunc(p1, angles), "r-", c='red')
#plt.xlim(1,2)
plt.show()

Optimized params: [ 3.61758709e-05  1.11037057e+02 -7.94918725e+01]

Оптимизация прошла на порядок лучше. Данных в статье не нашла(

In [30]:
#сделаем массив с шагом по величине торсионного угла CС
torses = np.arange(-180, 181, 12)
torses
Out[30]:
array([-180, -168, -156, -144, -132, -120, -108,  -96,  -84,  -72,  -60,
        -48,  -36,  -24,  -12,    0,   12,   24,   36,   48,   60,   72,
         84,   96,  108,  120,  132,  144,  156,  168,  180])
In [31]:
tors_energies = np.array([run_orca(tors.format(value))
                                 for value in torses])
In [32]:
tors_energies
Out[32]:
array([-79.57057865, -79.57008821, -79.56880146, -79.56719664,
       -79.56588321, -79.5653764 , -79.56588331, -79.56719512,
       -79.56880338, -79.57008835, -79.57057851, -79.57008821,
       -79.56880318, -79.56719666, -79.56588146, -79.5653764 ,
       -79.56588159, -79.56719683, -79.56880341, -79.57008835,
       -79.57057851, -79.57008821, -79.56880318, -79.56719666,
       -79.56588146, -79.5653764 , -79.56588159, -79.56719683,
       -79.56880341, -79.57008835, -79.57057851])
In [37]:
fitfunc = lambda p, x: p[0] * (1 + np.cos(p[1] * np.pi * np.radians(x) - p[2])) + p[3]  # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y  # Error function

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

#Plot it
plt.plot(torses, tors_energies, "ro", c='blue', alpha=0.5)
good_torses = np.linspace(-180, 180, 361)
plt.plot(good_torses, fitfunc(p1, good_torses), "r-", c='red')
#plt.xlim(1,2)
plt.show()

Optimized params: [ 2.59934980e-03  9.55058593e-01 -1.55516514e-05 -7.95705882e+01]

Минимумы находятся на отметках -180, -60, 60, 180 — с периодом в 120. Это описывает заторможенную конформацию этана.