Расчет энергии этана по теории функционала плотности с помощью GAMESS, перевод координат этана в формат pdb и подготовка pdb-файла:
gms ethane_zmt_dft.inp 1 >& ethane_zmt_dft.log
babel -igamout ethane_zmt_dft.log -opdb ethane.pdb
Ante_RED.pl ethane.pdb
В заголовке файла ethane-out.p2n указано, что этан имеет нулевой заряд и единичную мультиплетность, поэтому файл ethane-out.p2n был переименован в Mol_red1.p2n и были вычислены точечные заряды атомов этана:
mv ethane-out.p2n Mol_red1.p2n
RED-vIII.4.pl
@<TRIPOS>MOLECULE
LIG
8 7 1 0 1
SMALL
USER_CHARGES
@<TRIPOS>ATOM
1 C 0.763688 0.000000 -0.000000 C 1 LIG -0.0189 ****
2 H 1.156406 0.000000 -1.012074 H 1 LIG 0.0059 ****
3 H 1.156406 0.876482 0.506037 H 1 LIG 0.0059 ****
4 H 1.156406 -0.876482 0.506037 H 1 LIG 0.0059 ****
5 C -0.763688 0.000000 0.000000 C 1 LIG -0.0155 ****
6 H -1.156406 0.000000 1.012074 H 1 LIG 0.0056 ****
7 H -1.156406 0.876482 -0.506037 H 1 LIG 0.0056 ****
8 H -1.156406 -0.876482 -0.506037 H 1 LIG 0.0056 ****
@<TRIPOS>BOND
1 1 2 1
2 1 3 1
3 1 4 1
4 1 5 1
5 5 6 1
6 5 7 1
7 5 8 1
@<TRIPOS>SUBSTRUCTURE
1 LIG 1 **** 0 **** ****
Создадим описание этана в формате пакета программ GROMACS. При задании типов атомов и параметров функции Леннард-Джонса будем считать, что Ван-дер-Ваальсовые взаимодействия между атомами углерода разных молекул этана минимальны, т.к. атомы углерода почти полностью экранированы атомами водорода, что Ван-дер-Ваальсовый радиус водорода (sigma) известен, но не известен epsilon водорода. При описании молекулы укажем, что атомы через три связи не учитываются при расчете Ван-дер-Ваальсовых взаимодействий. При описании атомов укажем их точечные заряды, вычисленные ранее. При описании связей укажем их длины (в нанометрах), вычисленные ранее, и константы жесткости. Аналогично опишем углы и торсионные углы. Укажем пары атомов, которые не учитываются при расчете Ван-дер-Ваальсовых взаимодействий. При описании системы укажем координаты с 38 молекулами этана.
[ 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
ethane 3
[ 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
[ bonds ]
; ai aj funct b0 kb
1 2 1 0.1554 250000.0
1 3 1 0.1085 300000.0
1 4 1 0.1085 300000.0
1 5 1 0.1085 300000.0
2 6 1 0.1085 300000.0
2 7 1 0.1085 300000.0
2 8 1 0.1085 300000.0
[ angles ]
; ai aj ak funct phi0 kphi
;around c1
2 1 3 1 111.291 400.400
2 1 4 1 111.291 400.400
2 1 5 1 111.291 400.400
3 1 4 1 111.291 200.400
3 1 5 1 111.291 200.400
4 1 5 1 111.291 200.400
;around c2
1 2 6 1 111.291 400.400
1 2 7 1 111.291 400.400
1 2 8 1 111.291 400.400
6 2 7 1 111.291 200.400
6 2 8 1 111.291 200.400
7 2 8 1 111.291 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
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
[ pairs ]
; ai aj funct
3 6
3 7
3 8
4 6
4 7
4 8
5 6
5 7
5 8
[ System ]
first one
[ molecules ]
;Name count
ethane 38
Промоделируем испарение этана с помощью пакета программ GROMACS. Первая система соответствует газовой фазе, где расстояния между молекулами порядка 50 ангстрем. Вторая система имеет такую же плотность, как и жидкий этан. Создадим файлы, в которых epsilon водорода принимает различные значения, проведем моделирование динамики двух систем (для каждого файла) в программе grompp, вычислим энергию Ван-дер-Ваальсовых взаимодействий в каждой системе (для каждого файла) в программе g_energy, сохраним эти значения в файлы и определим разницу энергии Ван-дер-Ваальсовых взаимодействий между двумя системами (для каждого файла):
%%bash
echo Epsilon Delta_E_VdW
for i in {1..7}; do
ep=$(echo "scale=5; 1/$i/$i/$i" | bc -l )
sed "s/1.00000e-00/$ep/" ethane.top > v_${i}.top
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
echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10 vb_${i} -o eb_${i} > ethane_g_${i}.txt
echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10 v_${i} -o e_${i} > ethane_l_${i}.txt
gas=$(awk '$1=="LJ" {printf $3}' ethane_g_$i.txt)
liquid=$(awk '$1=="LJ" {printf $3}' ethane_l_$i.txt)
delta=$(echo $gas '-' $liquid | bc -l)
echo $ep $delta
done
При сравнении разницы энергии Ван-дер-Ваальсовых взаимодействий между двумя системами с экспериментальным значением энтальпии испарения этана (5.4 кДж/моль при 25°С) можно предположить, что значение epsilon водорода чуть больше, чем 0.037. Вкладом Кулоновских взаимодействий можно пренебречь.
Конвертируем соответствующие файлы в формат pdb в программе trjconv:
trjconv_d -f vb_3 -s vb_3 -o vb_3.pdb
trjconv_d -f v_3 -s v_3 -o v_3.pdb
Моделирование этана в жидкой фазе:
from IPython.display import HTML
HTML('<video controls><source src="v_3.mp4"></video>')
Моделирование этана в газовой фазе:
HTML('<video controls><source src="vb_3.mp4"></video>')