%%bash
cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
cat /usr/share/gromacs/top/oplsaa.ff/ffbonded.itp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
Получился следующий файл с топологией
(делала на kodomo, до того как он упал)
with open('./rem/topol.top','r') as fname:
for i in fname:
print(i.strip())
#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 3 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 1 2 7 1 1 2 8 1 6 2 7 1 6 2 8 1 7 2 8 1 [ dihedrals ] ; ai aj ak al funct 3 1 2 6 3 3 1 2 7 3 3 1 2 8 3 4 1 2 6 3 4 1 2 7 3 4 1 2 8 3 5 1 2 6 3 5 1 2 7 3 5 1 2 8 3 [ pairs ] ; list of atoms 1-4 ; 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
Загрузим для работы 5 конфигурационных файлов с разными алгоритмами контроля температуры:
- be.mdp - метод Берендсена для контроля температуры;
- vr.mdp - метод "Velocity rescale" для контроля температуры;
- nh.mdp - метод Нуза-Хувера для контроля температуры;
- an.mdp - метод Андерсена для контроля температуры;
- sd.mdp - метод стохастической молекулярной динамики.
%%bash
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp
--2021-05-24 13:12:48-- http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 93.180.63.127 Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:80... connected. HTTP request sent, awaiting response... 302 Found Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp [following] --2021-05-24 13:12:48-- https://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 1356 (1.3K) Saving to: ‘be.mdp’ 0K . 100% 114M=0s 2021-05-24 13:12:48 (114 MB/s) - ‘be.mdp’ saved [1356/1356] --2021-05-24 13:12:48-- http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 93.180.63.127 Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:80... connected. HTTP request sent, awaiting response... 302 Found Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp [following] --2021-05-24 13:12:48-- https://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 1427 (1.4K) Saving to: ‘vr.mdp’ 0K . 100% 102M=0s 2021-05-24 13:12:48 (102 MB/s) - ‘vr.mdp’ saved [1427/1427] --2021-05-24 13:12:48-- http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 93.180.63.127 Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:80... connected. HTTP request sent, awaiting response... 302 Found Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp [following] --2021-05-24 13:12:48-- https://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 1429 (1.4K) Saving to: ‘nh.mdp’ 0K . 100% 127M=0s 2021-05-24 13:12:48 (127 MB/s) - ‘nh.mdp’ saved [1429/1429] --2021-05-24 13:12:48-- http://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 93.180.63.127 Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:80... connected. HTTP request sent, awaiting response... 302 Found Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp [following] --2021-05-24 13:12:48-- https://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 1426 (1.4K) Saving to: ‘an.mdp’ 0K . 100% 128M=0s 2021-05-24 13:12:48 (128 MB/s) - ‘an.mdp’ saved [1426/1426] --2021-05-24 13:12:48-- http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 93.180.63.127 Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:80... connected. HTTP request sent, awaiting response... 302 Found Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp [following] --2021-05-24 13:12:48-- https://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 1441 (1.4K) Saving to: ‘sd.mdp’ 0K . 100% 133M=0s 2021-05-24 13:12:48 (133 MB/s) - ‘sd.mdp’ saved [1441/1441]
Сделаем входные файлы для молекулярно-динамического движка mdrun с помощью grompp
(ячейки ниже запускались на kodomo до того, как он упал (перенесла строки в другой ноутбук))
import os
sp = [i for i in os.listdir() if '.mdp' in i]
Ячейка ниже запускалась на kodomo до того, как он упал (перенесла строки в другой ноутбук):
for i in sp:
j = i.split('.')[0]+'.tpr'
k = i.split('.')[0]+'.trr'
m = i.split('.')[0]+'.pdb'
! grompp -f {i} -c etane.gro -p topol.top -o {j}
! mdrun -deffnm {i} -v -nt 1
! echo 0 | trjconv -f {k} -s {j} -o {m}
Посмотрим, что у нас получилось, в PyMol
from IPython.display import Image
Метод Берендсена
display(Image("./be.gif", format='png'))
Метод Андерсена
display(Image("./an.gif", format='png'))
Метод Нуза-Хувера
display(Image("./nh.gif", format='png'))
Velocity rescale
display(Image("./vr.gif", format='png'))
Метод стохастической MD
display(Image("./sd.gif", format='png'))
По траекториям мы видим следствие разного контроля температуры. Примечаательно, что при использовнии термостата Андерсена наша молекула практически "заморожена". В случае со стохастической MD результат противоположный - система перегревается. В случае с термостатом Берендсена молекула вращается неупорядоченно, не пребывая большей части времени в выгодной конформации, что тоже не вполне правильно.
import numpy as np
import matplotlib.pyplot as plt
def draw_plot(path, title):
mat = np.genfromtxt(path)
x, p, k = mat[:,0], mat[:,1], mat[:,2]
fig, ax = plt.subplots(figsize=(15,5))
ax.plot(x, p, label='потенциальная энергия', color='firebrick')
ax.plot(x, k, label='кинетическая энергия', color='cadetblue')
ax.set_title('Метод {}'.format(title), fontdict={'fontsize' : 18})
ax.set_xlabel('время (пс)', fontdict={'fontsize' : 14})
ax.set_ylabel('энергия (кДж/моль)', fontdict={'fontsize' : 14})
ax.legend()
draw_plot('./an.xvg', 'Андерсена')
Микровывод:
Как мы видим, никаких вращений метильных групп не происходила, только флуктации длины связи С-С (и углов). Термостат удерживает энергии на низких значениях - молекула действительно "подморожена"
draw_plot('./be.xvg', 'Берендсена')
Микровывод:
При использовании данного термостата молекула этана "не различала" свои состояния и метильные группы вращались равномерно. Мы видим, что флуктация энергии с течением времени симуляции уменьшается, приходя к относительно постоянному значению.
draw_plot('./vr.xvg', '"Velocity rescale"')
Микровывод:
Наблюдаются несколько более значимые флуктации энергии по сравнению с термостатами Андерсена и Берендсена. Тем не менее, "общий тренд" кажется постоянным.
draw_plot('./sd.xvg', 'Стохастической MD')
Микровывод:
Разброс значений энергий также, как и в Velocity Rescale побольше, чем в случае с термостатами Андерсена и Берендсена.
draw_plot('./nh.xvg', 'Нуза-Хувера')
Микровывод:
Наблюдается самый большой разброс значений энергии системы. Молекула исследует более широкий спектр энергетических состояний, периодически система будто бы остужается.
Рассмотрим распределение длинны связи С-С за время моделирования
def draw_hist(path, title):
mat = np.genfromtxt(path)
x, y = mat[:,0], mat[:,1]
fig, ax = plt.subplots(figsize=(8,4))
ax.bar(x, y, width = 0.0002, color='cadetblue')
ax.set_title('Метод {}'.format(title), fontdict={'fontsize' : 18})
ax.set_xlabel('длина связи (нм)', fontdict={'fontsize' : 14})
draw_hist('./an.1.xvg', 'Андерсена')
draw_hist('./be.1.xvg', 'Берендсена')
draw_hist('./vr.1.xvg', '"Velocity rescale"')
draw_hist('./sd.1.xvg', 'Стохастическая MD')
draw_hist('./nh.1.xvg', 'Нуза-Хувера')
В целом, все методы дали достаточно осмысленные результаты по длине связи С-С. По виду гистограммы, построенной для термостата Андерсена можно сказать, что он исследовал достаточно небольшое пространство длин связей. Также кажется, что распределение длин связей для метода "Velocity rescale" несколько сдвинуто в область более высоких значений.
import seaborn as sns
def draw_distribution(path, title):
mat = np.genfromtxt(path)
x, y = mat[:,0], mat[:,1]
fig, ax = plt.subplots(figsize=(8,4))
sns.histplot(y, color='cadetblue')
ax.set_title('Метод {}'.format(title), fontdict={'fontsize' : 18})
ax.set_xlabel('потенциальная энергия (кДж/моль)', fontdict={'fontsize' : 14})
draw_distribution('./an.xvg', 'Андерсена')
draw_distribution('./be.xvg', 'Берендсена')
draw_distribution('./vr.xvg', '"Velocity rescale"')
draw_distribution('./sd.xvg', 'Стохастическая MD')
draw_distribution('./nh.xvg', 'Нуза-Хувера')
Распределение Больцмана:
Из всех изображенных распределений наиболее приближенные картины дают те, что были получены с использованием методов Нуза-Хувера, Стохастической MD, "Velocity rescale".
Таблица зависимости быстродействия от алгоритма
import pandas as pd
# ns/day
dat = [5744.704, 5625.023, 5552.722, 4875.866, 5567.033]
pd.DataFrame(dat).rename({0 : 'an', 1 : 'be', 2 : 'vr', 3 : 'sd', 4 : 'nh'}).T
an | be | vr | sd | nh | |
---|---|---|---|---|---|
0 | 5744.704 | 5625.023 | 5552.722 | 4875.866 | 5567.033 |