Полезные ресурсы:
Уроки по работе с GROMACS находятся здесь.
Необходимые сведения о работе с GAMESS см. здесь.
Сведения о работе с Gnuplot см. здесь.
Введения о скриптовании в Bash здесь.
Ведения о awk здесь.
Цель задания: расчёт точечных зарядов на атомах этана. Построение файла топологии (этот файл содержит описание ковалентных и нековалентных взаимодействий). С помощью расчёта энтальпии испарения найти оптимальные параметры для описания Ван-дер-Ваальсовых взаимодействий.
Для начала определим точечные заряды. Для этого воспользуемся набором скриптов RED на perl. С помощью babel сделаем 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Если нет ошибок, то через какое-то время программа закончит работу. В директории Data-RED находится файл Mol_m1-o1.mol2 с координатами атомов и зарядами.
Теперь начнем создание файла описания молекулы в формате пакета программ GROMACS. Единица измерения расстояния в GROMACS - нанометр. Этот файл назовем et.top.
Как создавался файл:[ 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 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.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.
Здесь можно возразить: а как же 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
Создали описание молекулы с нуля. Чаще для этого используются программы с готовыми блоками.
Cледующая задача промоделировать испарение этана. Есть подготовленные состояния системы, первое соответствует газовой фазе, где расстояния между молекулами порядка 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
После этого необходимо провести для каждой системы молекулярную динамику с каждым файлом топологии. Дан файл с настройкам для динамики. Добавим в скрипт строчки для расчета:
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
Полученный скриптПолученные файлы для газа и для жидкой фазы. Для газа значения энергий VdW взаимодействий все порядка 10-4, тогда как у жидкости минимальный порядок 1. А кулоновские взаимодействия для обоих фаз намного меньше, чем VdW. Поэтому для оценки энтальпии испарения можно брать значения только для VdW взаимодействий в жидкой фазе. Epsblon водорода должна лежать в диапозоне от 0.016 до 0.037