Определим точечные заряды на атомах этана.
Для этого воспользуемся набором скриптов RED на perl.
С помощью babel сделаем pdb файл этана из результатов оптимизации из предыдущего практикума.
Добавим путь к скриптам в системный путь:
export PATH=${PATH}:/home/preps/golovin/progs/bin
Теперь с помощью скрипта Ante_RED.pl подготовим pdb файл:
Ante_RED.pl et.pdb
Переименуем p2n файл в Mol_red1.p2n и запустим RED:
RED-vIII.4.pl
По окончании работы программы находим файл Mol_m1-o1.mol2
с координатами атомов и зарядами.
Создадим файл с именем et.top описания молекулы в формате пакета программ GROMACS.
Единица измерения расстояния в GROMACS - нанометр.
Сначала задаём некоторые правила: [ defaults ].
Задаём типы атомов и параметры для функции Ленорда-Джонса - [ atomtypes ].
Будем считать, что в случае этана Ван-дер-Ваальсовое взаимодействие между атомами
углерода разных молекул минимально, так как углероды почти полностью экранированы
атомами водорода. Поэтому поставим для углерода некоторые параметры.
Ван-дер-Ваальсовый радиус водорода, т.е. сигма, нам известен - 120 pm
(webelements.com).
Итак у нас получилось, что в этом разделе всего одна переменная - это epsilon для водорода.
Дальше переходим непосредственно к описанию молекулы - [ moleculetype ].
Здесь описываем имя и указываем, что соседи через три связи не учитываются при
расчете Ван-дер-Ваальсовых взаимодействий. Это верно, так как мы включаем
это взаимодействие в торсионные углы.
Добавим атомы этана - [ atoms ].
Переходим к описанию связей - [ bonds ].
Константу жесткости и длину связи берём из предыдущего практикума.
Описание углов - то же самое. Силовые константы можно взять из примера:
[ 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
.........
Параметры торсионных углов возьмём из примера:
[ 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
.......
Создадим список пар атомов, которые не должны считаться при расчете VdW - [ pairs ].
Особенность расчета 1-4 взаимодействий подразумевает, что в профиле торсионного угла
участвует не только потенциал с cos, но и LJ отталкивание.
Для описания системы -[ System ]- воспользуемся любезно подготовленными г-ном
Головиным координатами с 38 молекулами этана:
[ System ] ; any text here first one [ molecules ] ;Name count et 38
В результате махинаций создали файл et.top.
Промоделируем испарение этана.
Воспользуемся двумя файлами с описаниями системы в газовой фазе с
расстояниями между молекулами порядка 50 ангстрем (box_big.gro)
и в жидкой фазе (box_38.gro).
Наша задача провести короткое моделирование динамики каждой из этих систем о определить
разницу в энергии VdW взаимодействий между системами, сравнить эту разницу с энтальпией
испарения этана. При Т=25 это значение равно 5.4 кДж/моль. Epsilon для водорода, кстати говоря, нам не известна.
Создадим 7 топологий с разными значениями epsilon. Будем использовать скрипт.
Теперь надо провести для каждой системы молекулярную динамику с каждым файлом топологии.
Для этого воспользуемся файлом с настройками для динамики.
Изменим скрипт для расчёта: et1.bash.
На основе полученных txt файлов установим среднее значение
энергии для каждого значения epsilon водорода:
для жидкости
для газа
Значения энергий VdW для газа на 4 порядка меньше, чем для жидкости, а кулоновские взаимодействия
значительно меньше VdW, поэтому при расчёте энтальпии испарения можно пренебречь VdW в
газовой фазе и кулоновскими взаимодействиями. На мой взгляд, для воспроизведения энтальпии
испарения этана epsilon должна находиться в пределах 0.00462 - 0.01562.