Молекулярная динамика биологических молекул в GROMACS


Цель данного занятия ознакомится с возможностями моделирования молекулярной динамики.

В этом занятии мы будем рользоваться пакетом молекулярной динамики Gromacs. Это программное обеспечение распостраняется под лицензией GPL, т.е. пользователь может скачать исходный код и свободен его изменять по своему усмотрению.

Общие положения

Типы файлов:

  • gro - файл с координатами системы.
  • top - файл с описанием ковалентных и нековалентных взаимодействий в молекулах.
  • mdp - файл с описанием параметров для работы молеклярно-механического движка.
  • tpr - файл для молеклярно-механического движка по сути есть объединение gro, top и mdp.
  • trr, xtc - файл с координатами после рассчёта.

Основные программы из пакета, которые будут использованны на занятии:

Программы запускаются в консоли, флаги запуска программ начинаются с -, например -f. Как правило после флага следует либо имя файла либо значение параметра. Смотрите примеры ниже.

  • editconf - манипуляция форматом координат и самими координатами. Пример: 

editconf -f my.gro -o my.pdb
  • genbox - наполнение ячейки растворителем.Пример: 

genbox -cp my.gro -cs mysolvent.gro -p my.top -o my_solvated.gro
  • genion - утилита для замены n молекул растворителя на ионы. 

genion -s my.tpr -np 10 -p my.top -o my_ions.gro 
-np это добавить 10 положительно заряженых ионов
  • grompp - объединение и проверка gro, top и mdp в tpr. 

grompp -f my.mdp -c my.gro -p my.top
  • mdrun - молеклярно-механический движок. На входе принимает tpr файл. 

mdrun -deffnm my
здесь параметр -deffnm означает, что выходные файлы будут называться как и входной файл, только с другими расщирениями

Настройка соединения с суперкомпьютером

Доступ к суперкомпьютеру возможен только по ssh ключам. Скопировать и настроить их использование: 

rm ~/.ssh
mkdir ~/.ssh
cp /home/preps/golovin/skif/* ~/.ssh
cat /home/preps/golovin/skif/config >> ~/.ssh/config
chmod 600 ~/.ssh/fbb-prac

Проверить соединение: 

ssh skif

Выбранный объект для практикума:

Порядок работы:

  1. Создайте рабочую директорию на удалённой машине (kodomo) с помощью WinSCP типа aleshin (чтобы удобнее было рабираться).
  2. Скачаем файлы:

    • дополнительной топологии для липида DPPC, dppc.itp.

    • параметры для липидов lipid.itp.

    • координаты одного липида dppc.gro.

    • Файл-заготовка тополгии системы b.top.

    • файл праметров для минимизации энергии em.mdp.

    • файл праметров для "утряски" воды pr.mdp pr.mdp.

    • файл праметров для молекулярной динамики md.mdp.

  1. в рабочей директории на kodomo на основе одного липида созадим ячейку с 64 липидами. 
    genconf -f dppc.gro -o b_64.gro -nbox 4 4 4

c помощью editconf (см. выше) создаем из b_64.gro файл b_64.pdb. Просмотрим результат в PyMol.

In [2]: 
#Для работы с Pymol: прокси и загрузка картинок:
from xmlrpclib import ServerProxy
from IPython.display import Image
import time

# launch pymol
import __main__
#__main__.pymol_argv = [ 'pymol', '-x' ]

__main__.pymol_argv = [ 'pymol', '-cp' ]  # to disable GUI
import pymol
from pymol import cmd
pymol.finish_launching()
In [3]: 
from IPython.display import Image

defaultImage = 'tmp\mypng1.png' 
def prepareImage(width=300, height=300, sleep=2, filename=defaultImage):
    ## To save the rendered image
    cmd.ray(width, height)
    cmd.png('tmp\mypng1.png')
    time.sleep(sleep)

def showPDBoneState(name):
    cmd.do('''
    reinit
    bg white
    set antialias, 2
    load ''' + name + '.pdb' + '''
    split_states ''' + name + '''
    hide everything, all
    select main, ''' + name + '_0001' + '''
    show spheres, main
    set sphere_scale, 0.2
    show sticks, main
    set stick_radius, 0.12
    zoom main
    ''')
In [6]: 
showPDBoneState('b_64')
prepareImage(); Image(defaultImage)
Out[6]:

Рис. 1. Созданный файл с 64 молекулами фосфолипида.

Как видно, все молекулы абсолютно идентичны и и идеально упорядочены.

  1. В текстовом редакторе в файле b.top установим правильное количество липидов в системе (64).
  2. Сделаем небольшой отступ в ячейке от липидов, чтобы добавить примерно 2500 молекул воды. 

editconf -f b_64.gro -o b_ec -d 0.5

  1. Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты молекул. 

grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v

выдача программы:

======= Making 3D domain decomposition 2 x 2 x 2

Steepest Descents:
Tolerance (Fmax) = 1.00000e+00
Number of steps = 1000000
Step= 0, Dmax= 2.0e-02 nm, Epot= 4.74007e+05 Fmax= 4.37970e+05, atom= 1842
Step= 1, Dmax= 2.0e-02 nm, Epot= 7.19293e+04 Fmax= 3.57933e+04, atom= 2588
Step= 2, Dmax= 2.4e-02 nm, Epot= 3.56795e+04 Fmax= 9.52568e+03, atom= 2993

...

Step= 50, Dmax= 1.6e-06 nm, Epot= 3.20591e+03 Fmax= 6.19379e+02, atom= 9155
Step= 51, Dmax= 1.9e-06 nm, Epot= 3.20911e+03 Fmax= 6.16887e+02, atom= 2765
Stepsize too small, or no change in energy.
Converged to machine precision,
but not to the requested precision Fmax < 1

Double precision normally gives you higher accuracy.

writing lowest energy coordinates.

Steepest Descents converged to machine precision in 52 steps,
but did not reach the requested Fmax < 1.
Potential Energy = 3.2059077e+03
Maximum force = 6.1937860e+02 on atom 915
Norm of force = 1.7839124e+02

Когда я выполнял это задание, не вчитался в комментарии, просто увидел причину ошибки (энергия выросла), поэтому не попробовал перезапустить программу по второму разу. Однако видно, что максимальная сила уменьшилась на три порядка.

  1. Добавим в ячейку молекулы воды типа spc. 

genbox -cp b_em -p b -cs spc216 -o b_s
  1. Проведём "утряску" воды: 

grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v

произошел взрыв системы, тогда: 

grompp -f em -c b_s -p b -o b_empr -maxwarn 1
mdrun -deffnm b_empr -v

и снова: 

grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v

  1. Переформатируем b_pr.gro и b_s.gro в pdb формат. И сравните визуально в PyMol изменения в системах:

In [9]: 
showPDBoneState('b_s')
prepareImage(); Image(defaultImage)
Out[9]:
In [10]: 
showPDBoneState('b_pr')
prepareImage(); Image(defaultImage)
Out[10]:

Рис. 2 & 3. Созданный файл с 64 молекулами фосфолипида и водой.

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

  1. Копируем файлы на суперкомпьютер: (разумеется, сперва надо настроить параметры (см. выше) 

cd ..
scp -r ./aleshin skif:_scratch/fbb
  1. Запускаем тестовое моделирование на суперкомпьютере. 

ssh skif
cd _scratch/fbb/aleshin
cp /home/users/golovin/progs/share/gromacs/top/residuetypes.dat .
grompp -f md -c b_pr -p b -o b_md -maxwarn 1
sbatch -n 4 -e error.log -o output.log -t 5 -p test impi /mnt/msu/users/fbbstudent/gmx-4.6.7-intel-impi-mkl-single/bin/mdrun_mpi -deffnm b_md -v

Запишите номер Вашей задачи. Просмотреть ход счёта можно в файле error.log 

error.log
Нажимая shift+. для перехода в конец файла. 
через пару минут процесс останавливается (лимит времени), ошибок не произошло.

Переходим дальше:

  1. Запускаем основное моделирование на суперкомпьютере. 

sbatch -N1 --ntasks-per-node=2 -e error-gpu.log -o output.log -t 350 -p gpu impi /opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi mdrun -testverlet -deffnm  b_md -v

Запишите номер Вашей задачи... К сожалению, Ломоносов перестал работать, задача ожидала в очереди, потому перезапуск файла error.log не произошел, если бы задача запустилась бы, снова можно было бы посмотреть 

less error.log
shift+.

и узнать время выполнения задачи. Также на Ломоносове есть возможность настроить отправку писем на почту о готовности работы, что замечательно, но надо разбираться.

в финале конвертируем файл ipython в html:

In [146]: 
%%bash
ipython nbconvert --to html pract8.ipynb
In []: