Получение координат атомов и их зарядов.
Для определения точечных зарядов воспользуемся набором скриптов RED на perl.С помощью babel был получен pdb файл этана из результатов оптимизации из предыдущего практикума:
В системный путь добавили путь к скриптам:
export PATH=${PATH}:/home/preps/golovin/progs/bin
babel -igamout et.log -opdb et.pdb
Теперь с помощью скрипта Ante_RED.pl подготовим pdb файл.
Ante_RED.pl et.pdb
Получились следующие файлы:
et.com, et-info.txt, et-out.p2n, et-out.pdb.
В заголовке в p2n файле указаны заряд и мультиплетность молекулы:
REMARK CHARGE-VALUE 0 (заряд молекулы = 0)
REMARK MULTIPLICITY-VALUE 1 (мультиплетность молекулы = 1)
Т.к. заряд и мультиплетность молекулы верны, переименуем p2n файл в Mol_red1.p2n.
Запустим RED:
RED-vIII.4.pl
Среди прочих получен файл Mol_m1-o1.mol2. Он содержит координаты атомов и заряды.
Создание файла топологии.
Следующий этап заключается в создании файла et.top с описанием молекулы в формате пакета программ GROMACS. Единица измерения расстояния в 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, константу жесткости для С-Н - из примера. et.log
!! Константа жесткости имеет разные размерности: gromacs kJ/(mol*nm), gamess eV/A.
Из википедии: В химии часто используется молярный эквивалент электронвольта. Если один моль электронов перенесён между точками с разностью потенциалов 1 В, он приобретает (или теряет) энергию Q = 96 485,3383(83) Дж, равную произведению 1 эВ на число Авогадро.
Пересчитаем значение k=0.563608:
0.563608*96 485,3383*10 = 543799,0855
[ bonds ]
; ai aj funct b0 kb
1 2 1 0.15298 543799.0855
1 3 1 0.10840 300000.0
1 4 1 0.10840 300000.0
1 5 1 0.10840 300000.0
2 6 1 0.10840 300000.0
2 7 1 0.10840 300000.0
2 8 1 0.10840 300000.0
Описание углов производим аналогично. Силовые константы были взяты из примера.
[ angles ]
; ai aj ak funct phi0 kphi
3 1 4 1 111.38 200.400
4 1 5 1 111.38 200.400
3 1 5 1 111.38 200.400
2 1 3 1 111.38 400.400
2 1 4 1 111.38 400.400
2 1 5 1 111.38 400.400
;around c2
1 2 6 1 111.38 400.400
6 2 8 1 111.38 200.400
6 2 7 1 111.38 200.400
7 2 8 1 111.38 200.400
1 2 7 1 111.38 400.400
1 2 8 1 111.38 400.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
3 1 2 8 1 0.0 0.62760 3
4 1 2 6 1 0.0 0.62760 3
4 1 2 7 1 0.0 0.62760 3
4 1 2 8 1 0.0 0.62760 3
5 1 2 6 1 0.0 0.62760 3
5 1 2 7 1 0.0 0.62760 3
5 1 2 8 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.
Моделирование испарения этана.
2 состояния системы: первое соответствует газовой фазе, где расстояния между молекулами порядка 50 ангстрем (box_big.gro), второе имеет такую же плотность, как и жидкий этан (box_38.gro).Необходимо провести короткое моделирование динамики каждой из этих систем и определить разницу в энергии VdW взаимодействий между системами, сравнив эту разницу с энтальпией испарения этана. При Т=25 это значение равно 5.4 кДж/моль.
Т.к. epsilon водорода неизвестна, молекулярная динамика и дальнейший расчет энергий были проведены для 7 топологий (по аналогии с занятием 4) с разными значениями epsilon с помощью скрипта.
Расчеты.
bash ./et_script.bashНа основе полученных txt файлов надо установить среднее значение энергии для каждого значения epsilon водорода и сравнить вклад кулоновских и VdW взаимодействий, а также оценить, в каком диапозоне должна лежать epsilon водорода, чтобы воспроизводилась энтальпия испарения этана.
Для молекул в жидкой фазе была получена следующая таблица:
epsilon average energy (kJ/mol)
LJ Coulomb
1.00000 -249.467 0.00259666
.12500 -28.8355 0.000163598
.03703 -8.25613 6.73105e-05
.01562 -2.62198 0.000228645
.00800 -2.91346 0.000261196
.00462 -2.26362 0.000166036
.00291 -0.980166 0.000129723
Для молекул в газовой фазе:
epsilon average energy (kJ/mol)
LJ Coulomb
1.00000 -0.000226728 3.87547e-05
.12500 -0.000158166 3.87606e-05
.03703 -0.000141376 3.89777e-05
.01562 -0.000134444 3.89762e-05
.00800 -0.000130805 3.89787e-05
.00462 -0.00012861 3.89867e-05
.00291 -0.000127178 3.89819e-05
Как видно из полученных данных, вклад кулоновских сил мал по сравнению с Ван-дер-Ваальсовыми.
Для газа значения энергий VdW взаимодействий все порядка 10-4,
тогда как у жидкости минимальный порядок 1.
Т.к. для обеих фаз кулоновские взаимодействия намного меньше, чем VdW,
то для оценки энтальпии испарения можно использовать значения только для VdW взаимодействий в жидкой фазе.
Чтобы воспроизводилась энтальпия испарения этана, epsilon водорода должна лежать в диапозоне от 0.01562 до 0.03703.