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

Полезные ресурсы:
Уроки по работе с GROMACS находятся здесь.
Необходимые сведения о работе с GAMESS см. здесь. 
Сведения о работе с Gnuplot см. здесь. 
Введения о скриптовании в Bash здесь. 
Ведения о awk здесь. 

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

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

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

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

   Ante_RED.pl et.pdb 

   Мультиплетность молекулы равна 1, заряд равен 0. Переименуем p2n-файл в Mol_red1.p2n и запустим RED:

   RED-vIII.4.pl

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

2. Теперь начнем создание файла описания молекулы в формате пакета программ GROMACS. Единица измерения расстояния в GROMACS - нанометр. Этот файл назовем et.top.

   Как создавался файл:
   В файлах этого типа комментарии находятся после ";". Итак первые две строчки, здесь задаются некоторые правила:

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

   [ bonds ]
   ;  ai    aj funct  b0       kb
        1   2   1  0.1554   150000.0 
        1   3   1  0.1085   180000.0
        1   4   1  0.1085   180000.0
        1   5   1  0.1085   180000.0
        2   6   1  0.1085   180000.0
        2   7   1  0.1085   180000.0
        2   8   1  0.1085   180000.0 

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

   [ angles ]
   ;  ai    aj    ak funct  phi0   kphi
   ;around c1
       3     1     4     1   109.500    200.400
       3     1     5     1   109.500    200.400
       4     1     5     1   109.500    200.400
       3     1     2     1   109.500    200.400
       4     1     2     1   109.500    200.400
       5     1     2     1   109.500    200.400
   ;around c2
       1     2     6     1   109.500    400.400
       1     2     7     1   109.500    400.400
       1     2     8     1   109.500    400.400
       6     2     7     1   109.500    200.400
       6     2     8     1   109.500    200.400
       7     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
       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

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

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

   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

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

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