запускаем 'ipython notebook' и 'pymol -R' из рабочей директории.
Форматирование задаётся Markdown: http://daringfireball.net/projects/markdown/basics
Работа проводится с помощью putty и обработки pymol. Команды для bash выделены табуляцией.
Изучение работы методов контроля температуры в GROMACS
#Для работы с 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()
запускаю ipython c Windows, в виду того, что порты на кодомо закрыты, это не предосудительно. У меня нет Bash, потому я просто копирую строки из скрипта ниже в putty, и они выполняются по мере готовности. Тем не менее, тэг "%%bash" я поставлю.
Изучим реализацию контроля температуры в молекулярной динамике на примере GROMACS. Объект исследования – это одна молекула этана.
Сначала подготовим файл координат и файл топологии. На основе .gro файла с 38 молекулами этана создадим индекс-файл, в котором будет группа из одной молекулы этана.
файл box_38.gro взят с http://kodomo.fbb.msu.ru/FBB/year_08/term6/box_38.gro.
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/box_38.gro%%bash
make_ndx -f box_38.gro -o 1.ndx
# Then type `r1` to choose the first residue
# Then type `q` to save and quit
Теперь создадим .gro файл с 1 молекулой и зададим ячейку.
При запуске ediconf нужно выбрать номер группы - выбираем одну молекулу - (3).
Затем зададим ячейку и расположим молекулу по центру ячейки.
%%bash
editconf -f box_38.gro -o et1.gro -n 1.ndx
# print 3
editconf -f et1.gro -o et.gro -d 2 -c
Построим файл топологии et.top для этана. При этом необходимо, в частности, подобрать типы атомов: /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
opls_135 12.01100 ; alkane CH3
opls_140 1.00800 ; alkane H.
ethane = """#include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
et 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 opls_135 1 ETH C1 1 -0.189 12.01
2 opls_135 1 ETH C2 2 -0.155 12.01
3 opls_140 1 ETH H1 3 0.0059 1.008
4 opls_140 1 ETH H2 4 0.0059 1.008
5 opls_140 1 ETH H3 5 0.0059 1.008
6 opls_140 1 ETH H4 6 0.0056 1.008
7 opls_140 1 ETH H5 7 0.0056 1.008
8 opls_140 1 ETH H6 8 0.0056 1.008
[ bonds ]
; ai aj funct b0 kb
1 2 1
1 4 1
1 5 1
2 6 1
2 7 1
2 8 1
[ angles ]
; ai aj ak funct phi0 kphi
;around c1
3 1 4 1
4 1 5 1
3 1 5 1
2 1 3 1
2 1 4 1
2 1 5 1
;around c2
1 2 6 1
6 2 8 1
6 2 7 1
7 2 8 1
1 2 7 1
1 2 8 1
[ dihedrals ]
; ai aj ak al funct
3 1 2 6 1
3 1 2 7 1
3 1 2 8 1
4 1 2 6 1
4 1 2 7 1
4 1 2 8 1
5 1 2 6 1
5 1 2 7 1
5 1 2 8 1
[ pairs ]
; ai aj funct
3 6
3 7
3 8
4 6
4 7
4 8
5 6
5 7
5 8
[ System ]
; any text here
first one
[ molecules ]
;Name count
et 1"""
with open('et.top', 'w') as w:
w.write(ethane)
так получили файл et.top (на рабочем компьютере), который был перекинут в рабочую директорию (putty) и использован sed:
sed -e 's/\r$//' et.top > et_sed.top но и так тестовый запуск (после скачивания http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp)
grompp -f be.mdp -c et.gro -p et.top -o et_test.tprили
grompp -f be.mdp -c et.gro -p et_sed.top -o et_test.tprне приводил к работе программы, которая ругалась на строки 44-52 файла et.top.
Лишь благодаря коллективному разуму (спасибо одному из немногих оставшихся однокурсникОВ) выход был найден с помощью замены параметра func с 1 на 3 в соответствующих строках ([ dihedrals ]).
%%bash
for i in $(echo be vr nh an sd); do
grompp -f $i.mdp -c et.gro -p et_sed.top -o et_$i.tpr
mdrun -deffnm et_$i -v -nt 1;
done
# sometimes you have to use ctrl+c
анализ результатов
Сперва посмотрим на pdb-файлы молекул.
%%bash
for i in $(echo be vr nh an sd); do
trjconv -f et_$i.trr -s et_$i.tpr -o et_$i.pdb;
done
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
''')
#showPDBoneState('et_be')
#showPDBoneState('et_vr')
showPDBoneState('et_nh')
#showPDBoneState('et_an')
#showPDBoneState('et_sd')
prepareImage(); Image(defaultImage)
Интересно наблюдать, как меняется набор состояний, проигрывая "запись" states в pymol для разных методов. Получилось, что число состояний в разных методах разное. Занятно, что в методе Нуза-Хувера изменяется лишь длина одной связи СН, а в остальных методах молекула скачет по ячейке, меняет углы и длины связей.
То есть можно ждать некой "замороженности" молекулы в методе Нуза-Хувера.
%%bash
for i in $(echo be vr nh an sd); do
g_energy -f et_$i.edr -o et_en_$i.xvg;
done
Так удалось посчитать потенциальную и кинетическую энергию, выбирая, значения 10 и 11 при выборе.
дальнейший анализ данных столкнулся с техническими трудностями (три блока команд далее) обойденными переходом на запуск ipython c kodomo.
%matplotlib inline
import matplotlib.pyplot as plt
import math
%lsmagic
%%bash
Если честно, уже порядком надоело, что не работает нормально magic. Я нашел похожее обсуждение на stackoverflow, установил Anaconda и обновил полностью ipython и т.д., прописал пути в PATH - безрезультатно. По какому принципу оказалось, что %lsmagic работает, а %matplotlib нет?..
###В итоге разобраться с этой проблемой дома у меня не получилось, поэтому я разобрался, как запускать ipython notebook на kodomo.
для себя прописываю на будущее ключевые команды:
перед запуском plink необходимо написать
set PATH=C:Files (x86)-0.63-RU-14;%PATH%
теперь мы можем использовать path из любой папки, например из рабочей директории, и далее следуем по инструкциям с http://kodomo.cmm.msu.ru/~golovin/ipynb/ipymol-kodomo.html
к сожалению запуск pymol с параметром -x сильно тормозит, но с параметром -cp все проходит гладко.
%matplotlib inline
import matplotlib.pyplot as plt
import math
УРА! все получилось! Я смогу использовать mathplotlib c kodomo и не заморачиваться с капризным Xserver!
%matplotlib inline
import matplotlib.pyplot as plt
import math
def parseXVG(filename): #будем обращаться к файлам по очереди, назовем каждый "filename" внутри этой функции
with open(filename, 'r') as f: # открываем на чтение созданный ранее файл
lines = filter(lambda x: not (x.startswith('#') or x.startswith('@')), f.readlines())
# с помощью лямбда-функции (чтобы не засорять код) фильтруем строки, пропуская те из них, которые начинаются с '#' или '@' все, что
# останется и есть интересующая нас информация.
[t, V, T] = zip(*[map(float, line.split()) for line in lines]) #здесь мы обработаем все строки из lines, разбив их
# на значимые числа и упаковав все в список: [время (пс), потенциальная энергия (кДж/моль), кинетическая энергия (кДж/моль)].
return [t, V, T]
def liftF(f, xs, ys):
return f(map(f, [xs, ys]))
def log10vec(vector):
return map(lambda x: math.log(x, 10), vector)
def showEnergy(filename, descr="", lpos="lower right"):
[x, y, z] = parseXVG(filename)
f = plt.figure(1)
f.text(0.5, 1, descr,
horizontalalignment='center', verticalalignment='top')
plt.subplot(141)
plt.title("Potential energy")
plt.scatter(x, y, c='b', alpha=0.7)
plt.xlabel("Time, ps")
plt.ylabel("Energy, kJ/mol")
plt.subplot(142)
plt.title("Kinetic energy")
plt.scatter(x, z, c='r', alpha=0.7)
plt.legend(loc=lpos)
plt.xlabel("Time, ps")
plt.subplot(143)
plt.title("log10 potential energy")
plt.scatter(x, log10vec(y), c='c', alpha=0.7)
plt.legend(loc=lpos)
plt.xlabel("Time, ps")
plt.ylabel("log10 energy, kJ/mol")
plt.subplot(144)
plt.title("log10 kinetic energy")
plt.scatter(x, log10vec(z), c='m', alpha=0.7)
plt.legend(loc=lpos)
plt.xlabel("Time, ps")
def showOneEnergy(listOfCoords, descr="", log=False, sparse=False, energy_h=25, energy_l=-.1): #получает список из списка времени(расположить по х) и энергий (по у).
labels = ('be', 'vr', 'nh', 'an', 'sd')
colours = ('r', 'g', 'b', 'y', 'c')
f = plt.figure(1)
f.text(0.5, 1, descr,
horizontalalignment='center', verticalalignment='top')
for i, pairedList in enumerate(listOfCoords): #работаем с списком из x (время) и y (энергия)
(x, y) = pairedList
plt.subplot(231+(i if i<4 else 5) if sparse else 151+i)
plt.title(labels[i])
plt.ylim(energy_l, energy_h) # тут пришлось покопаться и выбрать границы энергий для изображения всех графиков.
#Ясно, что оптимальным образом показать все не получится, поэтому ниже генерим целую серию рисунков
if log:
plt.scatter(x, log10vec(y), c=colours[i], alpha=0.7)
if i==0: plt.ylabel("log10 Energy")
else:
plt.scatter(x, y, c=colours[i], alpha=0.7)
if i==0: plt.ylabel("Energy, kJ/mol")
plt.xlabel("Time, ps")
plt.show()
filenames = ('et_en_be.xvg', 'et_en_vr.xvg', 'et_en_nh.xvg', 'et_en_an.xvg', 'et_en_sd.xvg')
coords = map(parseXVG, filenames)
getKinetic = lambda k: (k[0], k[2]) #обращается к списку и вытаскивает из него список с только первым и третьим элементами
getPotential = lambda k: (k[0], k[1]) #аналогично, но берет 1-й и 2-й элементы списка
Изобразим кинетические энергии, полученные разными методами. К сожалению, детальные различия показать не легко, потому придется генерировать серию графиков, соответствующих разным уровням энергии (все они в масштабе от нуля до i, где i от большего к меньшему: 5e17, 5e12, 5e10, 5e8, 5e6).
from matplotlib import pylab
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getKinetic, coords), descr="Kinetic energy, kJ/mol", energy_h=.5*1e18, energy_l=-.1*1e17)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getKinetic, coords), descr="Kinetic energy, kJ/mol", energy_h=.5*1e13, energy_l=-.1*1e12)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getKinetic, coords), descr="Kinetic energy, kJ/mol", energy_h=.5*1e11, energy_l=-.1*1e10)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getKinetic, coords), descr="Kinetic energy, kJ/mol", energy_h=.5*1e9, energy_l=-.1*1e8)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getKinetic, coords), descr="Kinetic energy, kJ/mol", energy_h=.5*1e7, energy_l=-.1*1e6)
Все выводы можно сделать уже из серии полученных изображений, но картина станет еще более ясной, если мы взглянем на логарифмический график энергий:
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getKinetic, coords), descr="Kinetic energy, kJ/mol", log=True)
аналогичные графики построим для потенциальной энергии, после чего сделаем выводы.
from matplotlib import pylab
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getPotential, coords), descr="Potential energy, kJ/mol", energy_h=.7*1e7, energy_l=-.1*1e6)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getPotential, coords), descr="Potential energy, kJ/mol", energy_h=.5*1e4, energy_l=-.1*1e3)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getPotential, coords), descr="Potential energy, kJ/mol", log=True, energy_h=8, energy_l=2)
Как же распределяется потенциальная энергия, какие значения чаще всего принимает?
def distPlotxxx(coords, type="scatter", cc_l=0, cc_h=3.6):
f = plt.figure(1)
f.text(0.5, 1, "Potential energy",
horizontalalignment='center', verticalalignment='top')
for i, pairedList in enumerate(coords):
(x, y) = pairedList
plt.subplot(151+i)
plt.title(labels[i])
if type=="hist":
plt.hist(y, color=colours[i], alpha=0.7)
else:
plt.scatter(x, y, color=colours[i], alpha=0.7)
plt.xlabel("Energy, kj/mol")
if i==0: plt.ylabel("Number of counts")
plt.ylim(cc_l, cc_h)
distCoords = map(parseXVG2d, filenames)
distPlotxxx(map(getPotential, coords), "hist", cc_l=0, cc_h=1000)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
Гистограммы здесь были призваны лишь показать форму распределения, что пригодится в конце страницы.
Ранее было сделано предсказание о "замороженности" молекулы этана при использованном методе Нуза-Хувера.
Теперь мы видим, что в данном случае кинетическая энергия практически постоянно равна приблизительно 1е8. На всякий случай проверил по файлу .xml, действительно, нигде ничего не перепутано, кинетическая энергия почти всегда более 100 млн кДж/моль. Насколько можно судить по графикам, это больше, чем медианное значение для других методов, но в отличие от Нуза-Хувера, значения кинетических энергий в них гораздо более разнообразны. Что касается потенциальной энергии молекулы в методе Нуза-Хувера, то она самая низкая среди рассмотренных методов и падает в самом начале измерений, после чего практически постоянна.
Вторым по аномальности можно считать метод Андерсена, где потенциальные энергии менее разнообразны, чем в трех оставшихся методах, а кинетическая энергия как бы скачет с уровня, аналогичного Нузу-Хуверу на более высокие уровни в три скачка, редко отклоняясь от принятого после перескока значения энергии.
К сожалению метод Берендсена считался очень недолго, после чего прекращал счет, потому мы имеем на порядок меньшее количество точек, но по полученным данным можно сказать, что в данном методе кинетические энергии молекулы в среднем меньше, чем в методе Velocity rescale и меньше, чем в методе стохастической молекулярной динамики.
Примерно то же можно сказать и о их потенциальной энергии, сделав оговорку, что в методе стохастической молекулярной динамики гораздо шире представлены и низкие и высокие энергии, распределенные во времени наиболее равномерно.
Теперь рассмотрим распределение длины связи С-С за время моделирования. Сначала создадим индекс файл с одной связью.
%%bash
echo -e "[ b ]\n1 2" > b.ndx
запустим утилиту по анализу связей g_bond.
%%bash
for i in $(echo be vr nh an sd); do
g_bond -f et_$i.trr -s et_$i.tpr -o bond_$i.xvg -d distance_$i.xvg -n b.ndx &> /dev/null;
done
запустим утилиту по анализу связей g_bond.
labels = ('be', 'vr', 'nh', 'an', 'sd')
filenames = map(lambda x: 'distance_%s.xvg' % x, labels)
colours = ('r', 'g', 'b', 'y', 'c')
def parseXVG2d(filename): # все аналогично прежней функции, только входной файл содержит не время и два типа энергии, а расстояние С-С.
with open(filename, 'r') as f:
lines = filter(lambda x: not (x.startswith('#') or x.startswith('@')), f.readlines())
[t, d] = zip(*[map(float, line.split()) for line in lines])
return [t, d] # [time (ps), distance (nm)]
def distPlot(coords, type="scatter", cc_l=0, cc_h=3.6):
f = plt.figure(1)
f.text(0.5, 1, "Distance",
horizontalalignment='center', verticalalignment='top')
for i, pairedList in enumerate(coords):
(x, y) = pairedList
plt.subplot(151+i)
plt.title(labels[i])
if type=="hist":
plt.hist(y, color=colours[i], alpha=0.7)
else:
plt.scatter(x, y, color=colours[i], alpha=0.7)
plt.xlabel("Distance, nm")
if i==0: plt.ylabel("Number of counts")
plt.ylim(cc_l, cc_h)
distCoords = map(parseXVG2d, filenames)
distPlot(distCoords, "hist", cc_l=0, cc_h=250)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
напоследок сменим масштаб оси, чтобы лучше изучить распределение.
distCoords = map(parseXVG2d, filenames)
distPlot(distCoords, "hist", cc_l=0, cc_h=35)
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
Что тут можно сказать. Странно, что связь СС совсем не колебалась в методе nh, странно, что длина связи в методе be изменяется раз в секунду.
Учитывая форму распределения Больцмана и распределения энергий в молекулах, посчитанных разными методами, очень трудно выбрать между методом vr и sd, но по совокупности признаков, как например, бОльшее разнообразие энергий и при этом бОльшая схожесть с распределением Больцмана, я бы предпочел указать метод стохастической молекулярной динамики, как способный наиболее реалистично поддерживать температуру в системе.
в финале конвертируем файл ipython в html:
ipython nbconvert --to html pract4.ipynb
%%bash
ipython nbconvert --to html pract7.ipynb