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

import subprocess
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

%matplotlib inline

Функция для запуска ORCA:

def run_orca(inp):
    with open('orca.inp', 'w') as outfile:
        outfile.write(inp)
    p = subprocess.Popen("/srv/databases/orca/orca orca.inp",  
                          shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out=p.communicate()[0]
    outlines = out.splitlines()
    for line in outlines:
        if "FINAL SINGLE" in line:
            energy = float(line.rstrip().split()[-1])
            return energy

Зависимость энергии молекулы этана от длины С-С связи

Имеем оптимизированную структуру этана в виде z-matrix. Создадим файлы для расчета энергии, варьируя длину одной из связей:

lengths = []
energies = []
step = 0.02
for n in range(-10, 11):
    length = 1.52986 + step*n
    lengths.append(length)
    inp = '''!HF RHF 6-31G
    * int 0 1
    C 0 0 0 0 0 0 
    C 1 0 0 %s 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
    *
    '''%(str(length))
    energies.append(run_orca(inp))

# x array
x_o = 1.52986 + step * np.arange(-10,11)
# y array
y_o = energies

Сохраним файл с энергиями:

with open('bond_en.txt', 'w') as f:
    f.write('value\tenergy\n')
    for v, en in zip(x_o, y_o):
        f.write('{}\t{}\n'.format(v, en))

Построим график:

#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,1, -79] # Initial guess for the parameters

def plot(x, y, start, end):
    p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
    print "Optimized params:", p1
    plt.plot(x_o, y_o, "ro", x_o, fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
    plt.xlim(start,end)
    plt.show()

plot(lengths, energies, 1.2, 1.8)

Optimized params: [  0.63991924   1.55699289 -79.08219401]

Зависимость энергии молекулы этана от валетного угла HCC

angles = []
angle_energies = []
step = 0.2
for n in range(21):
    angle = 109.2 + step*n
    angles.append(angle)
    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 %s 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
    *
    '''%(str(angle))
    angle_energies.append(run_orca(inp))

# x array
x_o = 109.2 + step * np.arange(21)
# y array
y_o = angle_energies

Сохраним файл с энергиями:

with open('angle_en.txt', 'w') as f:
    f.write('value\tenergy\n')
    for v, en in zip(x_o, y_o):
        f.write('{}\t{}\n'.format(v, en))

plot(angles, angle_energies, 108, 114)

Optimized params: [  4.10666297e-05   1.10890551e+02  -7.90815352e+01]

Зависимость энергии молекулы этана от торсионного угла CC

dihedrals = []
dih_energies = []
step = 12
for n in range(31):
    dihedral = -180 + step*n
    dihedrals.append(dihedral)
    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 %s
    H 2 1 6 1.08439 111.200 120
    H 2 1 6 1.08439 111.200 -120
    *
    '''%(str(dihedral))
    dih_energies.append(run_orca(inp))

# x array
x_o = -180 + step * np.arange(31)
# y array
y_o = dih_energies

Сохраним файл с энергиями:

with open('dihedral_en.txt', 'w') as f:
    f.write('value\tenergy\n')
    for v, en in zip(x_o, y_o):
        f.write('{}\t{}\n'.format(v, en))

fitfunc = lambda p, x: p[0] * (1 + np.cos(np.radians(3*x - p[1]))) + p[2]

plot(dihedrals, dih_energies, -190, 190)

Optimized params: [  2.29209805e-03   2.77352787e-03  -7.91975765e+01]

У изображенной выше функции 3 минимума, которые соответствуют трем заторможенным конформациям этана.

Обсуждение

Полученные параметры в целом похожи на те, что приведены в статье о разработке поля gaff (кроме жесткости связи). Отличия можно объяснить тем, что авторы использовали другие базисы для расчетов (MP2/6-31G*, MP4/6-311G**) и для получения конечных параметров использовали помимо ab inito вычислений также данные из кристаллических структур и из силового поля Amber.