В этом задании мы моделируем самосборку липидного слоя. Моделирование состоит из частей:
- Забрать с Kodomo разнообразные файлы топологий, координат и параметров
- Создать ячейку с N=64 липидами
- Преобразовать файлы ячеек в формат pdb с помощью editconf
- Скачать их и посмотреть в PyMol
- Оптимизировать систему, дабы сила контактов между молекулами была не экстремальной (?)
- Отметить максимальную силу контактов до и после оптимизации.
- Добавить в систему воды и «утрясти»
- Снова преобразовать в pdb и посмотреть на отличия
- Отправиться на Ломоносов и запустить там моделирование.
- Сначала тестовое, на 5 минут; если всё пошло хорошо — полноценное, на 350 минут.
Поехали.
Мы уже в правильной папке (Task9_mrkun). Грузим файлы.
%%bash
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/b.top
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/pr.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp
Создаём липидную ячейку и конвертируем её в pdb.
(Вывод GROMACS неинформативен, поэтому его мы опустим.)
%%bash
genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
editconf -f dppc.gro -o dppc.pdb
editconf -f b_64.gro -o b_64.pdb
Смотрим на файлы в PyMolnglview:
(А он умеет читать в том числе файлы .gro! Convenient.)
import nglview
import ipywidgets
from IPython.display import Image
view1 = nglview.show_structure_file('dppc.gro')
view1.clear_representations()
view1.add_representation(repr_type="licorice")
view1.center_view()
view1
#view1.download_image("dppc.png")
Image("dppc.png")
view2 = nglview.show_structure_file('b_64.gro')
view2.clear_representations()
view2.add_representation(repr_type="licorice")
view2.center_view()
view2
#view2.download_image("b_64.png")
Image("b_64.png")
Отредактировали b.top: там было указание на 34, а не 64 молекулы.
Отступаем от собственно липидов, чтобы добавлять воду не прямо в них, а на отдалении.
!editconf -f b_64.gro -o b_ec -d 0.5
Чиним плохие контакты:
!grompp -f em -c b_ec -p b -o b_em -maxwarn 2
!mdrun -deffnm b_em -v
Стартовая максимальная сила: Fmax = 4.37970e+05
Конечная максимальная сила: Fmax = 6.16887e+02
И это хорошо.
Добавляем воду, «трясём» её до базового состояния:
%%bash
genbox -cp b_em -p b -cs spc216 -o b_s
grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
Всё прошло хорошо, система не пошла вразнос.
Конвертируем b_pr.gro и b_s.gro в pdb.
!editconf -f b_pr.gro -o b_pr.pdb
!editconf -f b_s.gro -o b_s.pdb
Смотрим в nglview:
view3 = nglview.show_structure_file('b_s.gro')
view3.clear_representations()
view3.add_representation(repr_type="licorice")
view3.center_view()
view3
#view3.download_image("b_s.png")
Image("b_s.png")
view4 = nglview.show_structure_file('b_pr.gro')
view4.clear_representations()
view4.add_representation(repr_type="licorice")
view4.center_view()
view4
#view4.download_image("b_pr.png")
Image("b_pr.png")
Отправляемся на суперкомпьютер.
%%bash
cd ..
scp -r ./Task9_mrkun lom:_scratch/chem
Запускаем тестовое моделирование.
%%bash
ssh lom
cd _scratch/chem/Task9_mrkun
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 opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi \
mdrun -testverlet -deffnm b_md -v
Запускаем полное моделирование.
sbatch -N1 --ntasks-per-node=2 -e error_full.log -o output_full.log \
-t 350 -p gpu impi /opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi \
mdrun -testverlet -deffnm b_md -v
Номер полной задачи на gpu: 1643189
Анализ¶
!echo 2 | trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol
Посмотрим на образование бислоя глазами.
import pytraj
view5 = nglview.show_pytraj(pytraj.load("b_pbc_1.pdb"))
view5.clear_representations()
view5.add_representation(repr_type="licorice")
view5.center_view()
view5
Всё начинается с исходной регулярной структуры:
view5.render_image(frame=0)
view5._display_image()
Исходная структура быстро разрушается; уже к 20-30 кадру возникает гидрофобное ядро:
view5.render_image(frame=25)
view5._display_image()
К 50 кадру гидрофобное ядро отчётливо разделяется в бислой:
view5.render_image(frame=50)
view5._display_image()
Бислой сохраняется до конца моделирования, но около 70 кадра от него отделяется неспаренный фрагмент:
view5.render_image(frame=75)
view5._display_image()
view5.render_image(frame=100)
view5._display_image()
Посчитаем площадь, занимаемую одним липидом.
!echo 2 | g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline
import numpy as np
from __future__ import division
dims = np.loadtxt("box_1.xvg", skiprows=24)
dims[:6, :4]
plt.figure(figsize=(10,6))
# We want to compute per pair of lipids in a bilayer,
# so denominator is 32, not 64.
area = dims[:,2]*dims[:,3]/(4**3/2)
plt.plot(dims[:,0], area);
plt.show()
for stat in ("max", "mean", "median", "min", ):
print("{} area: {:f}".format(stat, getattr(np, stat)(area)))
Не сказать, что липиды укладываются стремительно, но зато укладываются хорошо: почти вдвое от стартового хаотичного состояния.
Посмотрим на площадь доступных гидрофильной и гидрофобной поверхностей.
!echo 2 2 | g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg -dt 100
hydro_areas = np.loadtxt("sas_b.xvg", skiprows=22)
hydro_areas[:6]
plt.figure(figsize=(10,6))
plt.plot(hydro_areas[:,0], hydro_areas[:,1],
color="red", label="Hydrophobic");
plt.plot(hydro_areas[:,0], hydro_areas[:,2],
color="blue", label="Hydrophilic");
plt.legend()
plt.show()
Площадь доступной гидрофобной поверхности стремительно снижается в первые ~5000 шагов, что соответствует увиденному выше образованию гидрофобного ядра. Это логично: упаковка гидрофобной поверхности в ядро существенно снижает энергию системы. В дальнейшем обе площади поверхностей продолжают снижаться, хотя и менее интенсивно: ядро обращается в бислой, в котором ещё меньше места для контактов с растворителем.
Считаем меру порядка бислоя в начале и в конце моделирования.
!g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
!g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X
order_start = np.loadtxt("ord_start.xvg", skiprows=12)
order_end = np.loadtxt("ord_end.xvg", skiprows=12)
order_end[:6]
plt.figure(figsize=(10,6))
plt.plot(order_start[:,0], order_start[:,3],
color="green", label="start");
plt.plot(order_end[:,0], order_end[:,3],
color="blue", label="end");
plt.legend()
plt.show()
Как и ожидалось, к концу моделирования билипидный слой существенно упорядоченнее, чем вскоре после начала.