Вычисление точечных зарядов и VdW параметров для молекулярной механики.
Суть задания состоит в расчёте точечных зарядов на атомах этана. Построение файла топологии, этот файл содержит описание ковалентных и нековалентных взаимодействий.
С помощью расчёта энтальпии испарения предлагается найти оптимальные параметры для описания Ван-дер-Ваальсовых взаимодействий.
-
Сначала необходимо определить точечные заряды. Для этого нужно воспользоываться набором скриптов RED на perl.
С помощью babel сделаем pdb-файл этана et.pdb из результатов оптимизации из предыдущего практикума. Добавим путь к скриптам в системный путь:
export PATH=${PATH}:/home/preps/golovin/progs/bin
Теперь с помощью скрипта Ante_RED.pl подготовим pdb файл:
Ante_RED.pl et.pdb
Мультиплетность молекулы равна 1, заряд равен 0.
Тогда переименуем p2n-файл в Mol_red1.p2n.
Теперь запустим RED:
RED-vIII.4.pl
Получили файл Mol_m1-o1.mol2 с координатами атомов и зарядами.
-
Начнем создание файла описания молекулы в формате пакета программ GROMACS. Единица измерения расстояния в GROMACS нанометр. Получили файл: et.top.
В первых двух строчках зададим некоторые правила:
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 0.5 0.8333
Дальше задаем типы атомов и собственно параметры для функции Леннорда-Джонса. Будем считать, что в случае этана Ван-дер-Ваальсовое взаимодействие между атомами углерода разных молекул минимально, так как углероды почти полностью экранированы атомами водорода. Поэтому поставим для углерода некоторые параметры.
Ван-дер-Ваальсовый радиус водорода (сигма) нам известен. Получилось, что в этом разделе имеем лишь одну переменную - 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 C 1 ETH C1 1 -0.0189 12.01
2 C 1 ETH C2 2 -0.0155 12.01
3 H 1 ETH H1 3 0.0059 1.008
4 H 1 ETH H2 4 0.0059 1.008
5 H 1 ETH H3 5 0.0059 1.008
6 H 1 ETH H4 6 0.0056 1.008
7 H 1 ETH H5 7 0.0056 1.008
8 H 1 ETH H6 8 0.0056 1.008
Переходим к описанию связей. Константу жесткости и длину связи возьмем из предыдущего занятия:
[ bonds ]
; ai aj funct b0 kb
1 2 1 0.1554 150000.0
1 3 1 0.1085 180000.0
1 4 1 0.1085 180000.0
1 5 1 0.1085 180000.0
2 6 1 0.1085 180000.0
2 7 1 0.1085 180000.0
2 8 1 0.1085 180000.0
Затем перейдем к описанию углов:
[ angles ]
; ai aj ak funct phi0 kphi
;around c1
3 1 4 1 109.500 200.400
3 1 5 1 109.500 200.400
4 1 5 1 109.500 200.400
3 1 2 1 109.500 200.400
4 1 2 1 109.500 200.400
5 1 2 1 109.500 200.400
;around c2
1 2 6 1 109.500 400.400
1 2 7 1 109.500 400.400
1 2 8 1 109.500 400.400
6 2 7 1 109.500 200.400
6 2 8 1 109.500 200.400
7 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
Теперь создадим список пар атомов, которые не должны считаться при расчете VdW. Особенность расчета 16-46 взаимодейс6твий подразумевает, что в профиле торсионного угла участвует не только потенциал с cos, но и LJ отталкивание. Это удобно для точной параметризации, но нам пока не надо.
Итак, добавляем список:
[ 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
В результате получаем файл et.top с полным описанием системы.
-
Мы создали описание молекулы с нуля. Наша следующая задача - промоделировать испарение этана. Воспользуемся двумя файлами с системами: первая соответствует газовой фазе, где расстояния между молекулами равны порядка 50 ангстрем, а вторая имеет такую же плотность, как и жидкий этан.
Наша задача заключается в том, чтобы провести короткое моделирование динамики каждой из этих систем, определить разницу в энергии VdW взаимодействий между системами и сравнить эту разницу с энтальпией испарения этана.
При Т=25 это значение равно 5.4 кДж/моль. Вспомним, что epsilon для водорода нам не известна.
По аналогии с предыдущим занятием создадим 7 топологий с разными значениями epsilon, после чего проведем для каждой системы молекулярную динамику с каждым файлом топологии, после чего посчитаем значения энергий (для чего воспользуемся утилитой g_energy).
Будем использовать скрипт:
#!/bin/bash
for i in {1..7}; do
###Создадим 7 топологий с разными значениями epsilon###
ep=$( echo "scale=5; 1/$i/$i/$i" | bc -l )
sed "s/1.00000/$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
###Посчитаем значения энергий с помощью утилиты g_energy###
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
<\pre>
После этого необходимо провести для каждой системы молекулярную динамику с каждым файлом топологии. Дан файл с настройкам для динамики. Добавим в скрипт строчки для расчета:
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
Теперь надо посчитать сами значения энергий, для этого воспользуемся утилитой 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
Скрипт был сохранен в файле.
-
В итоге мы получили txt файлы со средними значениями энергий для разных значений epsilon водорода.
Интересующие данные были сохранены в файлах epsflu.txt - для жидкой фазы и epsgas.txt - для газовой фазы.
Результаты показывают, что вклад VdW взаимодействий больше вклада кулоновских взаимодействий, особенно для системы жидкой фазы. Поэтому для оценки энтальпии испарения можно брать значения только для VdW взаимодействий в жидкой фазе. Epsblon водорода должна лежать в диапозоне от 0.016 до 0.037
<<Обратно на шестой семестр
<<Обратно на главную страницу
©Лелекова Мария,2011
|