|
Главная
I Семестр
II Семестр
III Семестр
IV Семестр
V Семестр
VI Семестр
Проекты
Обратная Связь
|
Вычисление точечных зарядов и VdW параметров для молекулярной механики
Начнем с того, что определим точечные заряды. Для этого воспользуемся набором скриптов
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 ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 0.5 0.8333
Далее заданы типы атомов и параметры для функции Ленорда-Джонса (считаем,
что в случае этана Ван-дер-Ваальсовое взаимодействие между атомами углерода
разных молекул минимально, так как углероды почти полностью экранированы атомами
водорода, поэтому поставим для углерода некоторые параметры).
Ван-дер-Ваальсовый радиус водорода, т.е. сигма, известен из многих источников,
см. webelements.com. В этом разделе всего одна переменная - это 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
Добавим атомы этана. В первом столбце идёт номер атома. На него мы будем ссылаться при
описании связей. Также необходимо добавить заряды для атомов, полученные с помощью RED.
[ 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
Переходим к описанию связей. Константу жесткости и длину связи C-C берём из
предыдущего занятия, а для связи С-H константу жёсткости берём так,
чтоб изменения значений для обеих связей были пропорциональны.
[ bonds ]
; ai aj funct b0 kb
1 2 1 0.1554 150000.0
1 3 1 0.1085 180000.0
.......
Описание углов, тоже самое. Силовые константы берём из примера.
[ angles ]
; ai aj ak funct phi0 kphi
;around c1
3 1 2 1 109.500 400.400
3 1 4 1 109.500 200.400
.........
;around c2
1 2 6 1 109.500 400.400
6 2 7 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 ]
; 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
Наша следующая задача промоделировать испарение этана. Даны два состояния
системы, первое соответствует газовой фазе, где расстояния между молекулами
порядка 50 ангстрем. Файл для газа.
Вторая система имеет такую же плотность как и жидкий этан.
Файл для жидкой фазы.
Наша задача провести короткое моделирование динамики каждой из этих
систем о определить разницу в энергии VdW взаимодействий между системами.
И сравнить эту разницу с энтальпией испарения этана. При Т=25 это
значение равно 5.4 кДж/моль. Вспомним, что epsilon для водорода нам
не известна. Создадим 7 топологий с разными значениями epsilon и
выполним для них расчёты.
Дан файл с настройкам для динамики.
Для поставленной задачи будем использовать bash скрипт:
make.bash.
for i in {1..7};do
#Cоздание 7 файлов топологий с
#разными значениями epsilon
ep=$( echo "scale=5; 1/$i/$i/$i" | bc -l )
sed "s/1.00000/$ep/" et.top > v_${i}.top
#Выполнение программы GROMACS, то есть проведение
#молекулярной динамики, для газовой и жидкой фаз
grompp_d -f md -c box_38 -p v_${i}.top -o v_${i} -maxwarn 1 && mdrun_d -deffnm v_${i} -v
grompp_d -f md -c box_big -p v_${i}.top -o vb_${i} -maxwarn 1 && mdrun_d -deffnm vb_${i} -v
#С помощью утилиты g_energy были
#посчитаны значения энергии
echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10 v_${i} -o e_${i} > v_${i}.txt
echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10 vb_${i} -o eb_${i} > vb_${i}.txt
#Запись значений epsilon и энергий
#VdW взаимодействий в один файл
echo -n "$ep " >> v_energy.txt
echo -n "$ep " >> vb_energy.txt
awk '/^LJ/{print $3}' v_${i}.txt >> v_energy.txt
awk '/^LJ/{print $3}' vb_${i}.txt >> vb_energy.txt
done
Можно конвертировать траекторию trr в pdb и посмотреть в PyMol.
trjconv_d -f v_1 -s v_1 -o v_1.pdb
Молекулы всех траекторий, кроме первой,
сразу же разлетались (иногда оставались летящие парами молекулы).
Полученные значения энергий VdW взаимодействий при соответствующих
значениях эпсилон находятся в файлах для газа и
для жидкости.
VdW взаимодействия в газообразном этане очень малы и имеют порядок не более 10-4,
в жидком этане энергия VdW имеет минимальный порядок 1.
Чтобы воспроизводилась энтальпия испарения этана epsilon водорода должна лежать в диапазоне от 0.01562 до 0.03703.
|