Kodomo

Пользователь

Учебная страница курса биоинформатики,
год поступления 2010

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

Отчёт по заданию должен появиться на сайте к следующему занятию.

Отчёт должен иметь ссылки на файлы с результатами счёта.

Традиционные ссылки на полезные ресурсы:

Вся работа по расчётам будет проходить на kodomo через терминал putty, а для работы с графическим выводом Gnuplot понадобится Xming.


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

  1. Начнем с того, что определим точечные заряды. Для этого воспользуемся набором скриптов RED на perl. С помощью babel сделаем pdb файл этана из результатов оптимизации из предыдущего практикума. Добавим путь к скриптам в системный путь: 

   1 export PATH=${PATH}:/home/preps/golovin/progs/bin

Теперь с помощью скрипта Ante_RED.pl подготовим pdb файл. 

   1 Ante_RED.pl et.pdb

Обратите внимание на заголовок в p2n файле, если Вас устраивает заряд и мультиплетность молекулы, то переименуйте Ваш p2n файл в Mol_red1.p2n. Запустите RED. 

   1 RED-vIII.4.pl

Обратите внимание на сообщения, если нет ошибок, то через какое-то время программа закончит работу и в директории Data-RED Вы сможете найти файл Mol_m1-o1.mol2 с координатами атомов и зарядами.

  1. Начнем создание файла описания молекулы в формате пакета программ GROMACS. Единица измерения расстояния в GROMACS нанометр. Пусть имя файла будет et.top. В дальнейшем я буду использовать это имя. В файлах этого типа комментарии находятся после ";". Итак первые две строчки, здесь мы задаём некоторые правила: 

   1 [ defaults ]
   2 ; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
   3 1               2               yes              0.5     0.8333

Дальше мы задаём типы атомов и собственно параметры для функции Ленорда-Джонса. Будем считать, что в случае этана Ван-дер-Ваальсовое взаимодействие между атомами углерода разных молекул минимально, так как углероды почти полностью экранированы атомами водорода. Поэтому поставим для углерода некоторые параметры. Ван-дер-Ваальсовый радиус водорода, т.е. сигма нам известен из многих источников, см. webelements.com. Итак у нас получилось, что в этом разделе всего одна переменная это epsilon для водорода. 

   1 [ atomtypes ]
   2 ; name      at.num  mass     charge ptype  sigma      epsilon
   3 H          1        1.008    0.0000  A   1.06908e-01  1.00000e-00
   4 C          6        12.01    0.0000  A   3.39967e-01  3.59824e-01

Дальше переходим непосредственно к описанию молекулы. Здесь мы описываем имя и указываем, что соседи через три связи не учитываются при расчете Ван-дер-Ваальсовых взаимодействий. Это верно так, как мы включаем это взаимодействие в торсионные углы: 

   1 [ moleculetype ]
   2 ; Name            nrexcl
   3 et            3

Добавим атомы этана. В первом столбце идёт номер атома. На него мы будем ссылаться при описании связей. Остальные столбца самоочевидные. Поправьте типы атомов и заряды. 

   1 [ atoms ]
   2 ;   nr  type  resnr  residue  atom   cgnr     charge       mass
   3      1   N      1    ETH      C1      1    -0.10      12.01
   4      2   N      1    ETH      C2      2    -0.10      12.01
   5      3   N      1    ETH      H1      3     0.0       1.008
   6      4   N      1    ETH      H2      4     0.0       1.008
   7      5   N      1    ETH      H3      5     0.0       1.008
   8      6   N      1    ETH      H4      6     0.0       1.008
   9      7   N      1    ETH      H5      7     0.0       1.008
  10      8   N      1    ETH      H6      8     0.0       1.008

Переходим к описанию связей. Константу жесткости и длину связи надо взять из занятия 4. 

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

дальше сами, связей должно быть 7

Описание углов, тоже самое. Силовые константы можно взять, как из примера. 

   1 [ angles ]
   2 ;  ai    aj    ak funct  phi0   kphi
   3 ;around c1
   4     3     1     4     1  109.500    200.400
   5 ..............
   6 ;around c2
   7     1     2     6     1   109.500    400.400
   8     6     2     8     1   109.500    200.400
   9 .........

должно быть 12 записей

Торсионные углы, возьмите параметры из примера: 

   1 [ dihedrals ]
   2 ;  ai    aj    ak    al funct  t0           kt      mult
   3     3    1     2     6      1  0.0      0.62760     3
   4     3    1     2     7      1  0.0      0.62760     3
   5 .......

должно быть 9 записей

Теперь давайте создадим список пар атомов которые не должны считаться при расчете VdW. Здесь можно возразить: а как же nrexcl=3 ? . Особенность расчета 1-4 взаимодействий подразумевает, что в профиле торсионного угла участвует не только потенциал с cos, но и LJ отталкивание. Это удобно для точной параметризации, но нам пока не надо. Итак добавляем список: 

   1 [ pairs ]
   2 ;  ai    aj funct
   3    3  6
   4    3  7
   5    3  8
   6    4  6
   7    4  7
   8    4  8
   9    5  6
  10    5  7
  11    5  8

Итак основное описание молекулы создано. Теперь переходим к описанию системы. Я подготовил для вас уже готовые координаты с 38 молекулами этана. Давайте укажем это в описании: 

   1 [ System ]
   2 ; any text here
   3 first one
   4 [ molecules ]
   5 ;Name count
   6  et    38

  1. Итак мы создали описание молекулы с нуля. Чаще для этого используются программы с готовыми блоками. Наша следующая задача промоделировать испарение этана. Тут всё просто, я подготовил два состояния системы, первое соответствует газовой фазе, где расстояния между молекулами порядка 50 ангстрем. Файл для газа. Вторая система имеет такую же плотность как и жидкий этан. Файл для жидкой фазы. Наша задача провести короткое моделирование динамики каждой из этих систем о определить разницу в энергии VdW взаимодействий между системами. И сравнить эту разнице с энтальпией испарения этана. При Т=25 это значение равно 5.4 кДж/моль. Вспомним, что epsilon для водорода нам не известна. Давайте по аналогии с занятием 4 создадим 7 топологий с разными значениями epsilon. Будем использовать скрипт: 

   1 #!/bin/bash
   2 for i in {1..7};do
   3     ep=$( echo "scale=5; 1/$i/$i/$i"  | bc -l )
   4     sed "s/тут надо что-то написать/$ep/" et.top > v_${i}.top
   5 done

Мы создали 7 файлов топологии. Теперь надо провести для каждой системы молекулярную динамику с каждым файлом топологии. Скачайте файл с настройкам для динамики. Добавим в скрипт строчки для расчета. 

   1  grompp_d -f md -c box_big -p v_${i}.top -o vb_${i} -maxwarn 1 && mdrun_d -deffnm  vb_${i} -v 
   2  grompp_d -f md -c box_38 -p v_${i}.top -o v_${i} -maxwarn 1 && mdrun_d -deffnm  v_${i} -v

Проведем расчет и комментируем строчки со счетом в скрипте. Желающие могут конвертировать траекторию trr в pdb и посмотреть в PyMol. 

   1 trjconv_d -f v_3 -s v_3 -o v_3.pdb

Теперь нам надо посчитать сами значения энергий, для этого воспользуемся утилитой g_energy. Эта утилита может работать в интерактивном режиме, но это не удобно в скрипте поэтому используем пере направление потока. Символ '\n' означает перенос строки: 

   1 echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10  vb_${i} -o eb_${i} > vb_${i}.txt
   2 echo -e "LJ-(SR)\nCoulomb-(SR)\n0" | g_energy -f -b 10  v_${i} -o e_${i} > v_${i}.txt
  1. На основе полученных txt файлов установите среднее значение энергии для каждого значения epsilon водорода. Сравните вклад кулоновских и VdW взаимодействий. Оцените в каком диапозоне должна лежать epsilon водорода, что бы воспроизводилась энтальпия испарения этана.
  2. * Установите точное значение epsilon для водорода.
  3. * Оцените вклад epsilon для углерода в разницу энергий между системами.