Начнем с того, что определим точечные заряды. Для этого воспользуемся набором скриптов RED на perl. С помощью babel сделаем pdb файл этана из результатов оптимизации из предыдущего практикума. Добавим путь к скриптам в системный путь:
export PATH=${PATH}:/home/preps/golovin/progs/bin
Теперь с помощью скрипта Ante_RED.pl подготовим pdb файл.
Ante_RED.pl et.pdbЗапустим 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
[ 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. Особенность расчета 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Основное описание молекулы создано. Теперь переходим к описанию системы:
[ System ] ; any text here first one [ molecules ] ;Name count et 38В результате получаем файл et.top с полным описанием системы.
Итак мы создали описание молекулы с нуля. Чаще для этого используются программы с готовыми блоками. Наша следующая задача промоделировать
испарение этана.
Для этого воспользуемся двумя файлами с системами, первая из которых соответствует газовой фазе, где расстояния между молекулами равны примерно 50 ангстрем (файл для газа), а вторая имеет такую же плотность, как и жидкий этан (файл для жидкой фазы). Кроме того, будем использовать файл с настройками для динамики. Проведем короткое моделирование динамики каждой из этих систем, определим разницу в энергии VdW взаимодействий между системами и сравним эту разницу с энтальпией испарения этана. При Т=25 это значение равно 5.4 кДж/моль. Вспомним, что epsilon для водорода нам не известна. По аналогии с предыдущим занятием создадим 7 топологий с разными значениями epsilon, после чего проведем для каждой системы молекулярную динамику с каждым файлом топологии, после чего посчитаем значения энергий (для чего воспользуемся утилитой g_energy). Скрипт для этой работы был сохранен в файле moldinscript.bash.
Полученные значения энергий VdW взаимодействий при соответствующих значениях эпсилон находятся в файлах для газа и для жидкости. VdW взаимодействия в газообразном этане очень малы и имеют порядок не более 10-4, в жидком этане энергия VdW имеет минимальный порядок 1. Чтобы воспроизводилась энтальпия испарения этана epsilop водорода должна лежать в диапазоне от 0.01562 до 0.03703.