Kodomo

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

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

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

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


Lifehack:

Запуск bash скриптов и просто программ из IPython Notebook :


  1. Вам предоставлена оптимизированная структура этана в виде z-matrix : 

   1  $DATA
   2 eth
   3 C1
   4  C   
   5  C      1   cc 
   6  H      2   ch   1   cchv 
   7  H      2   ch   1   cch   3   d1 0
   8  H      2   ch   1   cch   3   d2 0
   9  H      1   ch   2   cch   3   d3 0
  10  H      1   ch   2   cch   5   d3 0
  11  H      1   chv  2   cch   4   d3 0
  12 
  13 cc=1.52986
  14 ch=1.08439
  15 chv=1.08439
  16 cch=111.200
  17 cchv=111.200
  18 d1=120
  19 d2=-120
  20 d3=180
  21  $END

Как вы видите, вместо значений длин и углов связей стоят переменные. Наша цель состоит в том, что бы создать порядка 20 разных файлов для расчёта энергии в GAMESS с разными значениями по длине одной из связей. Для этого в представленных координатах одна связь или угол отличается названием переменной, но не её значением.

  1. Конечно можно сделать эти файлы вручную, но ожидается, что вы сделаете это средствами python . Можно попробовать автоматизировать процесс с помощью скрипта на bash. Bash это интерпретатор командной строки который автоматически запускается при присоединении к kodomo.

Составьте файл-заготовку для размножения, пусть имя файла будет et.inp. Для этого к координатам добавьте шапку для dft из предыдущего практикума. Только надо изменить информацию о типе входных координат: замените COORD=CART на COORD=ZMT. Проверьте работает ли Ваш файл-заготовка, т.е. запустите GAMESS если выходной файл не содержит ошибок, то переходим к следующему пункту.

  1. Теперь давайте создадим текстовый файл скрипта make_b.bash со следующим содержанием: 

   1 #!/bin/bash
   2 ### делаем цикл от -10 до 10 ##### 
   3 
   4 for i in {-10..10}; do 
   5 
   6 #### нам надо рассчитать новую длину связи #####
   7 #### с шагом 0.02 ангстрема,               #####
   8 #### воспользуемся калькулятором bc        #####
   9 #### и результат поместим в переменную nb  #####
  10 
  11      nb=$(echo "scale=5; 1.001 + $i/50" | bc -l)
  12 
  13 #### пролистаем файл et.inp и заменим указание переменной ###
  14 #### на новое значение и пере направим результат в файл   ###
  15 
  16      sed "s/cc=1.001/cc=$nb/" et.inp  > b_${i}.inp
  17 
  18 done

Поправьте скрипт так, что бы стартовая длина изменяемой связи соответствовала файлу et.inp. Запустите Ваш скрипт : 

   1  bash ./make_b.bash

Проверьте результат. У Вас должен быть 21 inp файл и в каждом разное значение для переменной сс.

  1. Давайте запустим расчёт для этих файлов, для этого перед строчкой с : 

   1 done

вставим запуск GAMESS: 

   1 gms b_${i}.inp 1 > b_${i}.log

Запустите скрипт и через какое-то время расчёт закончится. Теперь нам надо извлечь значение энергии из log файла. Удобно воспользоваться awk.
Сначала в нашем скрипте закомментируем запуск GAMESS поставив в начало строчки c gms # . Не пропустите это шаг! Добавим после за комментированной строчки вызов awk, при этом мы ищем строчку с TOTAL ENERGY и печатаем 4ое поле считая что поля разделены пробелами: 

   1 awk '/TOTAL ENERGY =/{print $4}'  b_${i}.log | tail -n 1

Запустите скрипт. На экране должно появится 21 значение энергии. Теперь удобно было бы выводить и значение длины связи, для этого добавьте перед вызовом awk распечатку переменой nb. Распечатаем переменную и несколько пробелов без переноса строки: 

   1 echo -n "$nb    "

Теперь запустите скрипт и если результат Вам подходит, т.е. есть две колонки цифр, то можно пере направить поток в файл : 

   1 bash ./make_b.bash > bond

  1. У вас есть зависимость энергии молекулы от длины одной связи. Эту зависимость можно построить в Excel. Вам предлагается сделать это в gnuplot или используйте matplotlib в IPython Notebook.

Для matplotlib: 

   1 %matplotlib inline
   2 import numpy as np
   3 import matplotlib.pyplot as plt
   4 from scipy import optimize


   1 a = np.loadtxt("bond")
   2 x_o=a[:,0]
   3 y_o=a[:,1]
   4 print "initial data:", y_o
   5 #function is  f(x)=k(b-x)^2 + a
   6 fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
   7 errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
   8 
   9 p0 = [1,1, -79] # Initial guess for the parameters
  10 p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
  11 print "Optimized params:", p1
  12 
  13 #Plot it
  14 plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
  15 plt.xlim(1.3,1.8)
  16 plt.show()

Для Gnuplot:

Запустите Xming->XLaunch. Выберите тип расположения окон, удобно использовать Multiple windows. Next. Выберите Start Program. Run Remote-> Putty. Дальше всё как обычно: kodomo, username.
Перейдите в рабочую директорию. Запустите Gnuplot: 

   1 gnuplot

Постройте зависимость энергии от длины связи, просто введите : plot "bond" У вас должен появиться график с точками похожими на параболу. Теперь нам надо найти коэффициенты в функции f(x)=a+k(x-b)^2 которые бы позволили наиболее близко описать наблюдаемую зависимость. Для этого воспользуемся возможностями Gnuplot. Сначала зададим функцию в развернутом виде, в строке gnuplot введём:  f(x)=a + k*x*x - 2*k*x*b + k*b*b  И зададим стартовые значения коэффициентов: 

   1 a=-80
   2 k=1
   3 b=1.5

Проведём подгонку коэффициентов под имеющиеся точки в файле bond:  fit f(x) "bond" via a,k,b  Сохраните значения коэффициентов. Постройте графики функции и значений энергии из GAMESS.  plot "bond", f(x) 

Сохраните изображение для отчёта. В отчёте обсудите причину не точного совпадения точек и функции.

  1. Проделайте аналогичные операции для валентного угла HCH, его значения должны изменяться от 109.2 до 113.2. Сохраните полученные коэффициенты в отчёт. Примечание: не перезаписывайте файл со значениями энергий они Вам нужны для отчёта. Сравните полученные константы с данными из статьи о GAFF. Укажите возможные причины расхождений ваших результатов и публикацией.
  2. Проделайте аналогичные операции для торсионного угла d3, его значения должны изменяться от -180 до 180 c шагом 12. Подгонку проводить не надо. Сохраните изображение и укажите в отчёте количество минимумов функции.
  3. * Увеличьте шаг до 0.1 ангстрема при расчёте связи. Постройте зависимость. Укажите какой функцией можно было бы аппроксимировать наблюдаемую зависимость.

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

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

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

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

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

   1 Ante_RED-1.5.pl et.pdb

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

   1 RED-vIII.51.pl

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

Найдите файл JOB2*log , а в файле раздел про "TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS". Сравните заряды из файла и полученные вами заряды с помощью RED (resp). Опишите разницу и предложите причину.

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

  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 для углерода в разницу энергий между системами.