Начнем с того, что определим точечные заряды. Для этого воспользуемся набором скриптов RED на perl. С помощью babel сделаем pdb файл этана из результатов оптимизации из предыдущего практикума. Добавим путь к скриптам в системный путь:
export PATH=${PATH}:/home/preps/golovin/progs/bin
Теперь с помощью скрипта Ante_RED.pl подготовим pdb файл.
Ante_RED.pl et.pdbОбратите внимание на заголовок в p2n файле, если Вас устраивает заряд и мультиплетность молекулы, то переименуйте Ваш 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.1525 250000.0
1 3 1 0.1085 300000.0
.......
дальше сами, связей должно быть 7
Описание углов, тоже самое. Силовые константы можно взять, как из примера.
[ 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
.........
должно быть 12 записей
Торсионные углы, возьмите параметры из примера:
[ 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
.......
должно быть 9 записей
Теперь давайте создадим список пар атомов которые не должны считаться при расчете 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
Итак мы создали описание молекулы с нуля. Чаще для этого используются программы с готовыми блоками. Наша следующая задача промоделировать испарение этана. Тут всё просто, я подготовил два состояния системы, первое соответствует газовой фазе, где расстояния между молекулами порядка 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
Мы создали 7 файлов топологии. Теперь надо провести для каждой системы молекулярную динамику с каждым файлом топологии. Скачайте
файл с настройкам для динамики. Добавим в скрипт строчки для расчета.
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
Проведем расчет и комментируем строчки со счетом в скрипте. trjconv_d -f v_3 -s v_3 -o v_3.pdbТеперь нам надо посчитать сами значения энергий, для этого воспользуемся утилитой 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
На основе полученных txt файлов установите среднее значение энергии для каждого значения epsilon водорода. Сравните вклад кулоновских и VdW взаимодействий. Оцените в каком диапозоне должна лежать epsilon водорода, что бы воспроизводилась энтальпия испарения этана.
* Установите точное значение epsilon для водорода.
* Оцените вклад epsilon для углерода в разницу энергий между системами.