Задание №8
Построим файл топологии для этана:
%%bash
cat ethane/ethane.top
Сгенерируем tpr файлы:
import subprocess
file_suffixes = ['be', 'vr', 'nh', 'an', 'sd']
method_names = ['Berendsen\'s Method', 'Velocity Rescale', 'Nose-Hoover MEthod', 'Andersen Method', 'Stochastic MD']
command = 'grompp -f {0}.mdp -c ethane.gro -p ethane.top -o ethane_{0}.tpr'
for suffix in file_suffixes:
subprocess.call(command.format(suffix), shell=True)
Для каждого из полученных файлов запустим mdrun и измерим время работы:
import time
command = 'mdrun -deffnm ethane_{0} -v -nt 1'
for method_name, suffix in zip(method_names, file_suffixes):
current_command = command.format(suffix)
start = time.time()
subprocess.call(current_command, shell=True)
print '{0}: {1} (seconds)'.format(method_names[ind], time.time() - start)
Конвертируем каждую из 5 систем в pdb:
command = 'echo 0 | trjconv -f ethane_{0}.trr -s ethane_{0}.tpr -o ethane_{0}.pdb'
for suffix in file_suffixes:
subprocess.call(command.format(suffix), shell=True)
Визуализируем получившиеся pdb:
from base64 import b64encode
from IPython.core.display import HTML, display
for method_name, suffix in zip(method_names, file_suffixes):
video = open('ethane_{0}.mp4'.format(suffix), 'rb').read()
video_encoded = b64encode(video)
html_code = '<h1>{0}</h1><video controls alt="PyMol Movie" src="data:video/mp4;base64,{1}" type="video/mp4">'.format(method_name, video_encoded)
display(HTML(html_code))
Видим, что при использовании метода Берендсена происходит активное вращение вокруг связи C-C. При velocity rescale и методы Нуза-Хувера также происходит вращение вокруг связи C-C, но менее интенсивное, чем при Берендсене. В целом velocity rescale похож на Нуза-Хувера. При методе Андерсена происходит лишь изменение длины связи C-C без вращений. И наконец стохастическая молекулярная динамик заставляет молекулу этана беспорядочно вращаться в пространстве. Кажется, что стохастическая молекулярная динамика показывает результаты имеющие минимальное отношение к действительности.
Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем:
command = 'echo -e \"10\n11\n\n\" | g_energy -f ethane_{0}.edr -o ethane_{0}_en.xvg -xvg none'
for suffix in file_suffixes:
subprocess.call(command.format(suffix), shell=True)
Построим график изменения энергий:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
for method_name, suffix in zip(method_names, file_suffixes):
fig = plt.figure(figsize=(15, 5))
ax = fig.add_subplot(111)
a = np.loadtxt('ethane_{0}_en.xvg'.format(suffix))
ax.plot(a[:,0], a[:,1], 'b', a[:,0], a[:,2], 'r')
ax.set_title(method_name)
blue_patch = mpatches.Patch(color='blue', label='Potential Energy')
red_patch = mpatches.Patch(color='red', label='Kinetic Energy')
ax.legend(handles=[blue_patch, red_patch])
plt.show()
Видим, что при методе Берендсена в отличие от других методов получается распределение энергии сильно похожее на распределение Больцмана.
Создадим индекс файл и запустим утилиту по анализу связей g_bond:
%%bash
echo -e "[ b ]\n1 2" > b.ndx
command = 'g_bond -f ethane_{0}.trr -s ethane_{0}.tpr -o bond_{0}.xvg -n b.ndx -xvg none'
for suffix in file_suffixes:
subprocess.call(command.format(suffix), shell=True)
Построим график распределения длин связей:
for method_name, suffix in zip(method_names, file_suffixes):
fig = plt.figure(figsize=(15, 5))
ax = fig.add_subplot(111)
a = np.loadtxt('bond_{0}.xvg'.format(suffix))
ax.set_title(method_name)
ax.bar(np.arange(len(a[:,0])), a[:,1])
plt.show()
Мы видим, что для методов Берендсена и Нуза-Хувера распределения длин связей больше всего похожи на нормальные. Помятуя также о том, что лишь при использовании метода Берендсена получается распределение энергий похожее на распределение Больцмана, можно сделать вывод, что в данной задаче для контроля температуры лучше всего подходит метод Берендсена, который к тому же еще и второй по быстродействию, уступаю методу Андерсена 0.05 секунды.