Занятие 5.

Вычисление точечных зарядов и VdW параметров для молекулярной механики
Суть задания состоит в расчёте точечных зарядов на атомах этана, построении файла топологии (этот файл содержит описание ковалентных и нековалентных взаимодействий). С помощью расчёта энтальпии испарения предлагается найти оптимальные параметры для описания Ван-дер-Ваальсовых взаимодействий.

Задание 1

Начнем с того, что определим точечные заряды. Для этого воспользуемся набором скриптов RED на perl.

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

export PATH=${PATH}:/home/preps/golovin/progs/bin

Теперь с помощью скрипта Ante_RED.pl подготовим pdb файл.

Ante_RED.pl et.pdb

Необходимо обратить внимание на заголовок в p2n файле, если нас устраивает заряд и мультиплетность молекулы (а нас это все усраивает:

REMARK CHARGE-VALUE 0 (заряд молекулы = 0)
REMARK MULTIPLICITY-VALUE 1 (мультиплетность молекулы = 1)),
то переименуем наш p2n файл в Mol_red1.p2n.

Запустим RED.

RED-vIII.4.pl

Программа не нашла ошибок и в директории Data-RED мы нашли файл Mol_m1-o1.mol2 с координатами атомов и зарядами.

Задание 2

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

Пусть имя файла будет et.top.

В файлах этого типа комментарии находятся после ";".

Итак, первые две строчки, в которых мы задаём некоторые правила:

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes              0.5     0.8333
Дальше мы задаём типы атомов и, собственно, параметры для функции Ленорда-Джонса.

Будем считать, что в случае этана Ван-дер-Ваальсовое взаимодействие между атомами углерода разных молекул минимально, так как углероды почти полностью экранированы атомами водорода.

Поэтому поставим для углерода некоторые параметры.

Ван-дер-Ваальсовый радиус водорода, т.е. сигма, нам известен из многих источников, см. webelements.com.

Итак, у нас получилось, что в этом разделе всего одна переменная - это epsilon для водорода.

[ 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   N      1    ETH      C1      1    -0.10      12.01
     2   N      1    ETH      C2      2    -0.10      12.01
     3   N      1    ETH      H1      3     0.0       1.008
     4   N      1    ETH      H2      4     0.0       1.008
     5   N      1    ETH      H3      5     0.0       1.008
     6   N      1    ETH      H4      6     0.0       1.008
     7   N      1    ETH      H5      7     0.0       1.008
     8   N      1    ETH      H6      8     0.0       1.008
Переходим к описанию связей. Константу жесткости и длину связи надо взять из занятия 4.

[ bonds ]
;  ai    aj funct  b0       kb
     1   2   1  0.1525   250000.0
     1   3   1  0.1085   300000.0
.......
дальше сами, связей должно быть 7
Описание углов производим аналогично. Силовые константы можно взять из примера.

[ angles ]
;  ai    aj    ak funct  phi0   kphi
;around c1
    3     1     4     1  109.500    200.400
..............
;around c2
    1     2     6     1   109.500    400.400
    6     2     8     1   109.500    200.400
.........
должно быть 12 записей
Торсионные углы, для них возьмем параметры из примера:

[ 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
.......
должно быть 9 записей
Теперь создадим список пар атомов, которые не должны считаться при расчете VdW.

(Здесь можно возразить: а как же nrexcl=3 ? .

Особенность расчета 1-4 взаимодействий подразумевает, что в профиле торсионного угла участвует не только потенциал с cos, но и LJ отталкивание. Это удобно для точной параметризации, но нам пока не надо.)

Итак, добавляем список:

[ pairs ]
;  ai    aj funct
   3  6
   3  7
   3  8
   4  6
   4  7
   4  8
   5  6
   5  7
   5  8
Основное описание молекулы создано. Ура-ура! Теперь переходим к описанию системы.

У нас есть уже готовые координаты с 38 молекулами этана.

Надо бы указать это в описании:

[ System ]
; any text here
first one
[ molecules ]
;Name count
 et    38
После строгого выполнения всего, о чем говорилось выше, получили вот какой файл et.top.

Задание 3

Мы создали описание молекулы с нуля.

Наша следующая задача - промоделировать испарение этана.

До этого было подготовлено два состояния системы (готовили, конечно, не мы, а преподаватель=)):

первое соответствует газовой фазе, где расстояния между молекулами порядка 50 ангстрем. Файл для газа.

Вторая система имеет такую же плотность, как и жидкий этан. Файл для жидкой фазы.

Наша задача состоит в том, чтобы провести короткое моделирование динамики каждой из этих систем и определить разницу в энергии VdW взаимодействий между системами, сравнив эту разницу с энтальпией испарения этана.

При Т=25 это значение равно 5.4 кДж/моль.

Вспомним, что epsilon для водорода нам не известна.

По аналогии с занятием 4 создадим 7 топологий с разными значениями epsilon.

Будем использовать скрипт:

#!/bin/bash
for i in {1..7};do
    ep=$( echo "scale=5; 1/$i/$i/$i"  | bc -l )
    sed "s/тут надо что-то написать/$ep/" et.top > v_${i}.top
done
Мы создали 7 файлов топологии. Теперь надо провести для каждой системы молекулярную динамику с каждым файлом топологии. Для этого мы скачали файл с настройкам для динамики.

Добавим в скрипт строчки для расчета: 

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
Проведем расчет и комментируем строчки со счетом в скрипте.

Далее конвертируем траекторию trr в pdb, при желании можно ее посмотреть в PyMol.

trjconv_d -f v_3 -s v_3 -o v_3.pdb

Полученный файл.

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

Эта утилита может работать в интерактивном режиме, но это не удобно, поэтому в скрипте используем перенаправление потока.

Символ '\n' означает перенос строки:

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
Снова выполнив все так, как было описано выше, получили конечный скрипт make_etan.bash.

Задание 4

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

Сравним вклад кулоновских и VdW взаимодействий.

Оценим, в каком диапозоне должна лежать epsilon водорода, чтобы воспроизводилась энтальпия испарения этана.

Получили следующую таблицу для молекул в жидкой фазе:


epsilon	average energy (kJ/mol)
             LJ	         Coulomb

epsilon	average energy (kJ/mol)
                  LJ	        Coulomb
1.00000	     -249.467	       0.00259666
.12500	     -28.8355	       0.000163598
.03703	     -8.25613	       6.73105e-05
.01562 	     -2.62198	       0.000228645
.00800 	     -2.91346	       0.000261196
.00462 	     -2.26362	       0.000166036
.00291 	     -0.980166	       0.000129723

Для молекул в газовой фазе:


epsilon	average energy (kJ/mol)
            LJ	               Coulomb
1.00000	-0.000226728	       3.87547e-05
.12500 	-0.000158166	       3.87606e-05
.03703 	-0.000141376	       3.89777e-05
.01562 	-0.000134444	       3.89762e-05
.00800 	-0.000130805	       3.89787e-05
.00462 	-0.00012861	       3.89867e-05
.00291 	-0.000127178	       3.89819e-05
Как видно из полученных данных, вклад кулоновских сил мал по сравнению с Ван-дер-Ваальсовыми. Для газа значения энергий VdW взаимодействий все порядка 10-4, тогда как у жидкости минимальный порядок 1.

Поскольку кулоновские взаимодействия для обоих фаз намного меньше, чем VdW, то для оценки энтальпии испарения можно брать значения только для VdW взаимодействий в жидкой фазе.

Чтобы воспроизводилась энтальпия испарения этана, epsilon водорода должна лежать в диапозоне от 0.01562 до 0.03703.