Учебный сайт

Бредихина Данилы

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

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

Суть задания состоит в определении констант ковалентных взаимодействий для молекулярной механики на основе квантово-химических расчётов.

Для оптимизированной структуры этана в виде z-матрицы создадим файлы для расчёта энергии в GAMESS, варьируя при этом длину одной из связей.

Для этого сделаем сначала файл-заготовку et.inp, состоящий из «шапки» для расчёта энергии с использованием теории функционала плотности (DFT) и структуры этана в виде z-матрицы.

In [14]: 
DFT_HEADER_ZMT=""" $CONTRL COORD=ZMT UNITS=ANGS   dfttyp=b3lyp RUNTYP=ENERGY $END
 $BASIS  GBASIS=N31 NGAUSS=6
 POLAR=POPN31 NDFUNC=1 $END
 $GUESS  GUESS=HUCKEL $END
 $system mwords=2 $end
"""
ETHANE_ZMT=""" $DATA
eth
C1
 C   
 C      1   cc 
 H      2   ch   1   cchv 
 H      2   ch   1   cch   3   d1 0
 H      2   ch   1   cch   3   d2 0
 H      1   ch   2   cch   3   d3 0
 H      1   ch   2   cch   5   d3 0
 H      1   chv  2   cch   4   d3 0

cc=1.52986
ch=1.08439
chv=1.08439
cch=111.200
cchv=111.200
d1=120
d2=-120
d3=180
 $END"""

with open('et.inp', 'w') as et:
    et.write(DFT_HEADER_ZMT)
    et.write(ETHANE_ZMT)

Теперь с помощью скрипта на bash создадим 21 файл. Длину связи в них, соответствующую переменной chv в z-матрице, будем менять с шагом 0.02 Å.

In [23]: 
%%bash
function change_var {
    # Usage: change_bond <input file> <variable to change> \
    #        <mean> <step> <n of iter / 2 - 1> <output prefix> <scale>
    for i in $(seq -$5 $5); do 
        nb=$(echo "scale=$7; $3 + $i*$4 / 1" | bc -l)
        sed -E "s/($2=)[0-9]+[.]?[0-9]*/\1$nb/" $1 > $6$i.inp;
    done
}
In [17]: 
%%bash
change_var et.inp cc 1.001 0.02 10 b_ 5

Запустим для созданных файлов программу GAMESS.

In [22]: 
GAMESS_SCRIPT="""#!/bin/bash
inp=$1
out=$(echo $inp | sed 's/\.inp/\.log/')
gms $inp 1 > $out
"""

with open('runGamess.sh', 'w') as rg:
    rg.write(GAMESS_SCRIPT)
In [24]: 
%%bash
chmod +x runGamess.sh
ls b_*.inp | xargs -n 1 ./runGamess.sh

Теперь извлечём значение энергии из log-файлов. С помощью программы awk будем выводить длину связи и значение энергнии для каждого из полученных файлов.

In [25]: 
%%bash
ls b_*.log | \
 xargs -n 1 awk -F '=' '/INPUT CARD>cc=/{printf $2} /TOTAL ENERGY =/{print $2}' | \
 sed -E 's/\s+/\t\t/' > bond.tsv

Изобразим графически завивисмоть энергии молекулы от длины одной связи.

In [5]: 
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

xs = []; ys = []
with open('bond.tsv', 'r') as bond:
    for line in bond:
        x, y = map(float, line.split())
        xs.append(x)
        ys.append(y)

plt.plot(xs, ys, 'k.')  # 'k.' means 'black dots'
plt.xlabel("C-C bond length, Angstrom")
plt.ylabel("Total energy, Ha")
Out[5]:
<matplotlib.text.Text at 0x7f6f3f277c50>

Попробуем описать полученные результаты в виде функции $ f(x)=ax^2+bx+c $.

In [6]: 
# Fitting
p = np.polyfit(xs, ys, 2)  # '2' is the degree of the fitting polynomial
xar = np.linspace(min(xs)-.05, max(xs)+.05, 100)

plt.clf()  # clear the figure
plt.plot(xs, ys, 'k.')
plt.plot(xar, np.polyval(p, xar), 'r-')  # A red solid line
plt.xlabel("C-C bond length, Angstrom")
plt.ylabel("Total energy, Ha")
Out[6]:
<matplotlib.text.Text at 0x7f6f3f962910>

Теперь аналогичные операции проделаем для одного из валентных углов HCC. Будем изменять его значения от 109,2° до 113,2°.

In [11]: 
%%bash
change_var et.inp cchv 111.2 0.2 10 a_ 5
In [12]: 
%%bash
ls a_*.inp | xargs -n 1 ./runGamess.sh
ls a_*.log | \
 xargs -n 1 awk -F '=' '/INPUT CARD>cchv=/{printf $2} /TOTAL ENERGY =/{print $2}' | \
 sed -E 's/\s+/\t\t/' > angle.tsv

Изобразим полученные результаты графически.

In [13]: 
xs = []; ys = []
with open('angle.tsv', 'r') as bond:
    for line in bond:
        x, y = map(float, line.split())
        xs.append(x)
        ys.append(y)

plt.plot(xs, ys, 'k.')  # 'k.' means 'black dots'

p = np.polyfit(xs, ys, 2)
xar = np.linspace(min(xs)-.05, max(xs)+.05, 100)
plt.plot(xar, np.polyval(p, xar), 'b-')

plt.xlabel(r"HCC angle, $ ^{\circ} $")
plt.ylabel("Total energy, Ha")

Теперь аналогичные манипуляции проделаем для одного из торсионных углов. Будем изменять значения угла d3 от -180° до 180° с шагом 12.

In [14]: 
%%bash
change_var et.inp d3 0 12 15 t_ 0
In [15]: 
%%bash
ls t_*.inp | xargs -n 1 ./runGamess.sh
ls t_*.log | \
 xargs -n 1 awk -F '=' '/INPUT CARD>d3=/{printf $2} /TOTAL ENERGY =/{print $2}' | \
 sed -E 's/\s+/\t\t/' > torsion.tsv

Изобразим полученные результаты на графике.

In [17]: 
xs = []; ys = []
with open('torsion.tsv', 'r') as bond:
    for line in bond:
        x, y = map(float, line.split())
        xs.append(x)
        ys.append(y)

plt.plot(xs, ys, 'k.')

plt.xlabel(r"d3 torsion angle, $ ^{\circ} $")
plt.ylabel("Total energy, Ha")
Out[17]:
<matplotlib.text.Text at 0x7f6f3f962390>

На графике отчётливо различимы минимумы функции.

Теперь увеличим шаг до 0.1 Å при расчёте длины связи и построим зависимость энергии от значения длины связи.

In [19]: 
%%bash
change_var et.inp cc 1.001 0.1 9 b2_ 5
In [20]: 
%%bash
ls b2_*.inp | xargs -n 1 ./runGamess.sh
ls b2_*.log | \
 xargs -n 1 awk -F '=' '/INPUT CARD>cc=/{printf $2} /TOTAL ENERGY =/{print $2}' | \
 sed -E 's/\s+/\t\t/' > bond2.tsv
In [21]: 
xs = []; ys = []
with open('bond2.tsv', 'r') as bond:
    for line in bond:
        x, y = map(float, line.split())
        xs.append(x)
        ys.append(y)

plt.plot(xs, ys, 'k.')  # 'k.' means 'black dots'

plt.xlabel("C-C bond length, Angstrom")
plt.ylabel("Total energy, Ha")

Наблюдаемую зависимость можно попытаться аппроксимировать гиперболой.

(Построим график $ y $.)

In [29]: 
xs = []; ys = []
with open('bond2.tsv', 'r') as bond:
    for line in bond:
        x, y = map(float, line.split())
        xs.append(x)
        ys.append(y)

newX = map(lambda x: 1/x, xs)

plt.plot(newX, ys, 'k.')
p = np.polyfit(newX, ys, 3)
xar = np.linspace(min(newX)-.05, max(newX)+.05, 100)
plt.plot(xar, np.polyval(p, xar), 'b-')

plt.xlabel(r"Inverse of C-C bond length, $Angstrom^{-1}$")
plt.ylabel(r"Total energy, $Ha$")
plt.show()

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

Суть этой части задания состоит в расчёте точечных зарядов на атомах этана. С помощью расчёта энтальпии испарения найдём оптимальные параметры для описания Ван-дер-Ваальсовых взаимодействий.

Определим точечные заряды. Для этого воспользуемся набором скриптов RED на Perl. С помощью babel сделаем .pdb файл этана из результатов оптимизации.

In []: 
%%bash
sed 's/ENERGY/OPTIMIZE/' et.inp > et2.inp
gms et2.inp 1 >& et2.log

babel -igamout et2.log -opdb et.pdb

Подготовим .pdb файл для дальнейших расчётов.

In []: 
%%bash
Ante_RED.pl et.pdb
mv et-out.p2n Mol_red1.p2n
RED-vIII.4.pl
ls Data-RED/Mol_m1-o1.mol2

Создадим файл описания молекулы в формате пакета программ GROMACS.

In [2]: 
gromacs_data = """[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes              0.5     0.8333
[ atomtypes ]
; name      at.num  mass     charge ptype  sigma      epsilon
H          1        1.008    0.0000  A   1.06908e-01  1.00000e-00
C          6        12.01    0.0000  A   3.39967e-01  3.59824e-01
[ moleculetype ]
; Name            nrexcl
et            3
[ atoms ]
;   nr  type  resnr  residue  atom   cgnr     charge       mass
     1   C      1    ETH      C1      1     0.0      12.01
     2   C      1    ETH      C2      2     0.0      12.01
     3   H      1    ETH      H1      3     0.0       1.008
     4   H      1    ETH      H2      4     0.0       1.008
     5   H      1    ETH      H3      5     0.0       1.008
     6   H      1    ETH      H4      6     0.0       1.008
     7   H      1    ETH      H5      7     0.0       1.008
     8   H      1    ETH      H6      8     0.0       1.008
[ bonds ]
;  ai    aj funct  b0       kb
     1   2   1  0.1525   250000.0
     1   3   1  0.1085   300000.0
     1   4   1  0.1085   300000.0
     1   5   1  0.1085   300000.0
     2   6   1  0.1085   300000.0
     2   7   1  0.1085   300000.0
     2   8   1  0.1085   300000.0
[ angles ]
;  ai    aj    ak funct  phi0   kphi
;around c1
    3     1     4     1  109.500    200.400
    4     1     5     1  109.500    200.400
    3     1     5     1  109.500    200.400
    2     1     3     1  109.500    200.400
    2     1     4     1  109.500    200.400
    2     1     5     1  109.500    200.400
;around c2
    1     2     6     1  109.500    200.400 
    6     2     8     1  109.500    200.400  
    6     2     7     1  109.500    200.400 
    7     2     8     1  109.500    200.400 
    1     2     7     1  109.500    200.400 
    1     2     8     1  109.500    200.400 
[ dihedrals ]
;  ai    aj    ak    al funct  t0           kt      mult
    3    1     2     6      1  0.0      0.62760     3
    3    1     2     7      1  0.0      0.62760     3
    3    1     2     8      1  0.0      0.62760     3
    4    1     2     6      1  0.0      0.62760     3
    4    1     2     7      1  0.0      0.62760     3
    4    1     2     8      1  0.0      0.62760     3
    5    1     2     6      1  0.0      0.62760     3
    5    1     2     7      1  0.0      0.62760     3
    5    1     2     8      1  0.0      0.62760     3
[ pairs ]
;  ai    aj funct
   3  6
   3  7
   3  8
   4  6
   4  7
   4  8
   5  6
   5  7
   5  8
[ System ]
; any text here
first one
[ molecules ]
;Name count
 et    38"""

with open('et.top', 'w') as w:
    w.write(gromacs_data)

Итак, мы создали описание молекулы этана.

Теперь промоделируем испарение этана. Воспользуемся для этого двумя файлами, первый из которых соответствует состоянию системы в газовой фазе, а второй – жидкому этану. Проведём короткое моделирование динамики каждой из этих систем и определим разницу в энергии VdW взаимодействий между системами.

Так как значение epsilon для водорода нам не известно, создадим 7 топологий с разными значениями epsilon. Затем проведём для каждой системы молекулярную динамику с каждым файлом топологии. (Для этого потребуется также файл с настройкам для динамики.)

In [4]: 
%%bash
for i in {1..7}; do
    ep=$( echo "scale=5; 1/$i/$i/$i"  | bc -l )
    sed "s/1.00000e-00/$ep/" et.top > v_${i}.top
    grompp_d -f md -c box_big -p v_${i}.top -o vb_${i} -maxwarn 1 && mdrun_d -deffnm  vb_${i} -v 
    grompp_d -f md -c box_38 -p v_${i}.top -o v_${i} -maxwarn 1 && mdrun_d -deffnm  v_${i} -v;
done

Посчитаем сами значения энергий, воспользовавшись для этого утилитой g_energy.

In []: 
%%bash
for i in {1..7}; do
    echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10  vb_${i} -o eb_${i} > vb_${i}.txt
    echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10  v_${i} -o e_${i} > v_${i}.txt;
done

На основе полученных .txt файлов установим среднее значение энергии для каждого значения epsilon водорода.

Кулоновские взаимодействия вносят нулевой вклад (рассматриваемые молекулы не заряжены). Поэтому расчёт произведём для взаимодействий Леннарда-Джонса.

In [5]: 
%%bash
for i in {1..7}; do
    echo -ne $i"\t"
    echo $(awk '$1=="LJ" {print $3}' vb_$i.txt) "-" $(awk '$1=="LJ" {print $3}' v_$i.txt) | bc -l;
done > evaporation.tsv

cat evaporation.tsv

1	267.227773529
2	21.468842159
3	2.432578814
4	3.914165734
5	3.752159367
6	.263324558
7	3.206632987

Сравним полученные значения с энтальпией испарения этана. При \( Т=25^{\circ} \) это значение равно \( 5.4 кДж/моль \). Наиболее близкая к нему величина была получена при \( \epsilon = 4 \).

Ссылки на файлы

  1. et.inp
  2. bond.tsv
  3. angle.tsv
  4. torsion.tsv
  5. bond2.tsv
  6. et.top
  7. Отчёт в iPython Notebook.
△