1. Начнем с того, что определим точечные заряды. Для этого воспользуемся набором скриптов RED на perl. С помощью babel сделаем pdb файл этана из результатов оптимизации из предыдущего практикума:
babel -ilog et.log -opdb et.pdb
Добавим путь к скриптам в системный путь:
export PATH=${PATH}:/home/preps/golovin/progs/bin
Теперь с помощью скрипта Ante_RED.pl подготовим pdb файл:
Ante_RED.pl et.pdb
Нас устраивает заряд и мультиплетность молекулы:
REMARK CHARGE-VALUE 0
REMARK MULTIPLICITY-VALUE 1
Поэтому переменовываем p2n файл в Mol_red1.p2n и запускаем RED:
RED-vIII.4.pl
В результате в директории Data-RED появится файл Mol_m1-o1.mol2 с координатами атомов и зарядами.
2. Начнем создание файла описания молекулы в формате пакета программ GROMACS.
В первых двух строчках мы задаём некоторые правила:
[ 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
Добавим атомы этана. В первом столбце идёт номер атома. На него мы будем ссылаться при описании связей.
[ 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. (7 связей)
[ bonds ]
; ai aj funct b0 kb
1 2 1 0.1525 250000.0
1 3 1 0.1085 300000.0
.......
Описание углов производим аналогично. Силовые константы можно взять из примера. (12 записей)
[ 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 ........
Торсионные углы, для них возьмем параметры из примера: (9 записей)
[ 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.
Здесь можно возразить: а как же 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
По окончанию заполнения получаем файл et.top .
3. Наша следующая задача промоделировать испарение этана.
Даны два состояния системы, первое соответствует газовой фазе, где расстояния между молекулами порядка 50 ангстрем.Файл для газа.
Вторая система имеет такую же плотность как и жидкий этан. Файл для жидкой фазы.
Наша задача провести короткое моделирование динамики каждой из этих систем о определить разницу в энергии VdW взаимодействий между системами.
И сравнить эту разницу с энтальпией испарения этана. При Т=25 это значение равно 5.4 кДж/моль.
Вспомним, что epsilon для водорода нам не известна. Поэтому необходимо создать 7 топологий с разными значениями epsilon и выполнить для них расчёты.
Будем использовать скрипт, в основе которого будет лежать следующее:
Мы создали 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
Проведем расчет и прокомментируем строчки со счетом в скрипте. Скрипт.
4. На основе полученных txt файлов необходимо установить среднее значение энергии для каждого значения epsilon водорода.
Для газа значения энергий VdW взаимодействий все порядка 10-4, тогда как у жидкости минимальный порядок 1. Для оценки энтальпии испарения можно брать значения только для VdW взаимодействий в жидкой фазе, пренебрегая кулоновскими взаимодейтсвиями, потому что для обоих фаз они намного меньше, чем VdW. Полученные файлы: для газа, для жидкости.
Оцените в каком диапозоне должна лежать epsilon водорода, что бы воспроизводилась энтальпия испарения этана:
Epsilon водорода должна лежать в диапозоне от 0.016 до 0.037 .