Начнем с того, что определим точечные заряды. Для этого воспользуемся набором скриптов 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
В директории Data-RED мы можем найти файл Mol_m1-o1.mol2 с координатами атомов и зарядами.Создадим файл описания молекулы в формате пакета программ GROMACS. Единица измерения расстояния в GROMACS нанометр. Пусть имя файла будет et.top. В файлах этого типа комментарии находятся после ";".
В первых двух строчках мы задаём некоторые правила [ defaults ].
Дальше мы задаём типы атомов и собственно параметры для функции Ленорда-Джонса [ atomtypes ]. Будем считать, что в случае этана Ван-дер-Ваальсовое взаимодействие между атомами углерода разных молекул минимально, так как углероды почти полностью экранированы атомами водорода. Поэтому поставим для углерода некоторые параметры. Ван-дер-Ваальсовый радиус водорода, т.е. сигма нам известен из многих источников, см. webelements.com. Итак у нас получилось, что в этом разделе всего одна переменная это epsilon для водорода.
Дальше переходим непосредственно к описанию молекулы [ moleculetype ]. Здесь мы описываем имя и указываем, что соседи через три связи не учитываются при расчете Ван-дер-Ваальсовых взаимодействий. Это верно так, как мы включаем это взаимодействие в торсионные углы.
Добавим атомы этана [ atoms ]. В первом столбце идёт номер атома. На него мы будем ссылаться при описании связей.
Переходим к описанию связей. Константу жесткости и длину связи C-C берём из предыдущего занятия, а для связи С-H константу жёсткости берём так, чтоб изменения значений для обеих связей были пропорциональны.
[ bonds ] ; ai aj funct b0 kb 1 2 1 0.1525 250000.0 1 3 1 0.1085 300000.0 .......Описание углов, тоже самое. Силовые константы берём из примера.
[ 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 ]. Здесь можно возразить: а как же nrexcl=3 ? . Особенность расчета 1-4 взаимодействий подразумевает, что в профиле торсионного угла участвует не только потенциал с cos, но и LJ отталкивание. Это удобно для точной параметризации, но нам пока не надо.Итак основное описание молекулы создано. Теперь переходим к описанию системы. Я подготовил для вас уже готовые координаты с 38 молекулами этана. Давайте укажем это в описании:
[ System ] ; any text here first one [ molecules ] ;Name count et 38
Готовый файл et.top.
Наша следующая задача промоделировать испарение этана. Даны два состояния системы, первое соответствует газовой фазе, где расстояния между молекулами порядка 50 ангстрем. Файл для газа. Вторая система имеет такую же плотность как и жидкий этан. Файл для жидкой фазы. Наша задача провести короткое моделирование динамики каждой из этих систем о определить разницу в энергии VdW взаимодействий между системами. И сравнить эту разницу с энтальпией испарения этана. При Т=25 это значение равно 5.4 кДж/моль. Вспомним, что epsilon для водорода нам не известна. Создадим 7 топологий с разными значениями epsilon и выполним для них расчёты.
Дан файл с настройкам для динамикиБудем использовать скрипт:
#!/bin/bash ### Создадим отдельные папки для ### ### хранения файлов, относящихся ### ### к топологии, жидкой и газовой ### фазам ### mkdir v-top mkdir v mkdir vb for i in {1..7}; do ### Расчёт значения epsilon ### ep=$( echo "scale=5; 1/$i/$i/$i" | bc -l ) ### Замена значения epsilon и запись ### ### результата в новый файл ### sed "s/1.00000/$ep/" et.top > v-top/v_${i}.top ### Выполнение программы GROMACS, ### ### то есть проведение молекулярной ### ### динамики, для газовой и жидкой фаз ### grompp_d -f md -c box_big -p v-top/v_${i}.top -o vb/vb_${i} -maxwarn 1 && mdrun_d -deffnm vb/vb_${i} -v grompp_d -f md -c box_38 -p v-top/v_${i}.top -o v/v_${i} -maxwarn 1 && mdrun_d -deffnm v/v_${i} -v ### Выполнение утилиты g_energy для ### ### расчёта энергий VdW и Кулоновских ### ### взаимодействий ### echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10 vb/vb_${i} -o vb/eb_${i} > vb/vb_${i}.txt echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10 v/v_${i} -o v/e_${i} > v/v_${i}.txt ### Запись значений epsilon и энергий ### ### VdW взаимодействий в один файл ### echo -n "$ep " >> vb/vb_energ.txt echo -n "$ep " >> v/v_energ.txt awk '/^LJ/{print $3}' vb/vb_${i}.txt >> vb/vb_energ.txt awk '/^LJ/{print $3}' v/v_${i}.txt >> v/v_energ.txt doneМожно конвертировать траекторию trr в pdb и посмотреть в PyMol.
trjconv_d -f v_3 -s v_3 -o v_3.pdb
Интересной мне показалась траектория v_1.pdb.На основе полученных txt файлов необходимо установить среднее значение энергии для каждого значения epsilon водорода.
Для газа значения энергий VdW взаимодействий все порядка 10-4, тогда как у жидкости минимальный порядок 1 (для жидкости, для газа). А кулоновские взаимодействия для обоих фаз намного меньше, чем VdW. Поэтому для оценки энтальпии испарения можно брать значения только для VdW взаимодействий в жидкой фазе.
Чтобы воспроизводилась энтальпия испарения этана, epsblon водорода должна лежать в диапозоне от 0.01562 до 0.03703.