Aznauryan's page

Суть задания состоит в расчете точечных зарядов на атомах этана. Построение файла топологии, этот файл содержит описание ковалентных и нековалентных взаимодействий. С помощью расчета энтальпии испарения предлагается найти оптимальные параметры для описания Ван-дер-Ваальсовых взаимодействий.

Определим точечные заряды. Для этого воспользуемся набором скриптов 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

Дальше мы задаем типы атомов и собственно параметры для функции Ленорда-Джонса. Будем считать, что в случае этана Ван-дер-Ваальсовое взаимодействие между атомами углерода разных молекул минимально, так как углероды почти полностью экранированы атомами водорода. Поэтому поставим для углерода некоторые параметры. Ван-дер-Ваальсовый радиус водорода, т.е. сигма нам известен из многих источников. Итак у нас получилось, что в этом разделе всего одна переменная это 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.

[ bonds ]
;  ai    aj funct  b0       kb
     1   2   1  0.1525   250000.0
     1   3   1  0.1085   300000.0
.......

Описание углов, тоже самое.

[ 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
.........

Торсионные углы:

[ 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

В итоге получаем файл et.top.

Итак мы создали описание молекулы с нуля. Следующая задача промоделировать испарение этана. Даны два состояния системы, первое соответствует газовой фазе, где расстояния между молекулами порядка 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

Теперь нам надо посчитать сами значения энергий, для этого воспользуемся утилитой 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

Полученный скрипт.

На основе полученных файлов попробуем установить среднее значение энергии для каждого значения epsilon водорода. Для газа:

1.00000    -0.00022681
0.12500     -0.000158047
0.03703     -0.00014137
0.01562     -0.000134441
0.00800     -0.000130803
0.00462     -0.000128609
0.00291     -0.000127177

Для жидкой фазы:

1.00000    -248.051
0.12500     -27.8758
0.03703     -6.17972
0.01562     -2.15295
0.00800     -0.980635
0.00462     -1.3682
0.00291     -2.13419

Таким образом, для газа значения энергий VdW взаимодействий все порядка 10-4, тогда как у жидкости минимальный порядок 1. А кулоновские взаимодействия для обоих фаз намного меньше, чем VdW. Поэтому для оценки энтальпии испарения можно брать значения только для VdW взаимодействий в жидкой фазе. Чтобы воспроизводилась энтальпия испарения этана epsilon водорода должна лежать в диапазоне от 0.01562 до 0.03703.

Вычисление точечных зарядов и VdW параметровдля молекулярной механики

`Вычисление точечных зарядов и VdW параметров для молекулярной механики`