Изучение работы методов контроля температуры в молекулярной динамике¶
%%bash
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
Файл et.top был дополнен связями и углами, типами атомов из ffbonded.itp. Тип функционала для торсионных углов - 3, а не 1.
!cat et.top
Даны 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
!ls
Из-за необходимости ручного ввода группы для преобразования в pdb, пришлось запускать команды ниже вручную.
types = ['be', 'vr', 'nh', 'an', 'sd']
for t in types:
mdp = t + '.mdp'
tpr = t + '.tpr'
trr = t + '.trr'
pdb = t + '.pdb'
com1 = f"grompp -f {mdp} -c etane.gro -p et.top -o {tpr}"
com2 = f"mdrun -deffnm {t} -v -nt 1"
com3 = f"trjconv -f {trr} -s {tpr} -o {pdb}"
print(com1)
print(com2)
print(com3)
print("")
!ls
from IPython.display import Image
Посмотрим на полученные результаты в PyMol. Они были визуализированы локально, сохранены покадрово, а затем собраны с помощью ffmpeg.
Метод Берендсена
display(Image("be.gif", format='png'))
Метод "Velocity rescale"
display(Image("vr.gif", format='png'))
Метод Нуза-Хувера
display(Image("nh.gif", format='png'))
Метод Андерсена
display(Image("an.gif", format='png'))
Метод стохастической молекулярной динамики
display(Image("sd.gif", format='png'))
Вывод: наименьшая подвижность наблюдается для контроля температуры методом Андерсена, наибольшая - для стохастической молекулярной динамики. Для метода "Velocity rescale" создается впечатление, что противолежащие водороды вращаются вместе, не заслоняя друг друга.
Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем.¶
Снова запустим команды ниже вручную, потому необходимо дополнительно вводить параметры выдачи: 10 - потенциальная энергия и 11 - кинетическая.
for t in types:
com1 = f"g_energy -f {t}.edr -o {t}.xvg -xvg none"
print(com1)
import numpy as np
import matplotlib.pyplot as plt
def plot_energy(short, long):
data_in = np.genfromtxt(short+'.xvg')
x, pot, kin = data_in[:,0], data_in[:,1], data_in[:,2]
fig, ax = plt.subplots(figsize=(10,5))
ax.plot(x, pot, label='Потенциальная E')
ax.plot(x, kin, label='Кинетическая E')
ax.set_title(f'Метод {long}')
ax.set_xlabel('Время (пс)')
ax.set_ylabel('Энергия (кДж/моль)')
ax.legend()
methods = ['Берендсена', '"Velocity rescale"', 'Нуза-Хувера',
'Андерсена', 'стохастической молекулярной динамики']
for t, m in zip(types, methods):
plot_energy(t, m)
Энергия этана при методе Андерсена оказалось наименьшей, что объясняет низкую подвижность молекулы при визуализации. Флуктуация энергии в методе Берендсена со временем уменьшается, а величина энергии остается более-менее постоянной. В остальных случаях энергии меняются в больших диапазонах. Для метода Нуза-Хувера наблюдается значительный разброс: энрегия оказывается как очень маленькой, так и очень большой. Однако чаще энергия находится на уровне менее 75 кДж/моль, а пики больше 150 хоть и встречаются, но не очень часто.
Распределение длины C-C связи¶
Сначала создадим индекс файл с одной связью.
!echo '[ bonds ]' > b.ndx
!echo '1 2' >> b.ndx
for t in types:
com1 = f"g_bond -f {t}.trr -s {t}.tpr -o bond_{t}.xvg -n b.ndx -xvg none"
!$com1
!ls | grep bond
def plot_bonds(short, long):
data_in = np.genfromtxt('bond_'+short+'.xvg')
x, l = data_in[:,0], data_in[:,1]
fig, ax = plt.subplots(figsize=(10,5))
ax.bar(x, l, width=0.0005)
ax.set_title(f'Метод {long}')
ax.set_xlabel('Длина связи, нм')
methods = ['Берендсена', '"Velocity rescale"', 'Нуза-Хувера',
'Андерсена', 'стохастической молекулярной динамики']
for t, m in zip(types, methods):
plot_bonds(t, m)
Для метода Андерсена диапазон колебаний длин связи оказался наименьшим, что согласуется с наименьшей энергией и общей заторможенностью структуры. Также заметно, что распределение для данного метода напоминает главное здание МГУ, то есть имеет 2 заметных пика, кроме центрального. У метода Берендсена наблюдается узкий пик для длины связи, что объясняется узким диапазоном энергий. У остальных методов распределения по связям сопоставимы, хотя у метода Нуза-Хувера пик всё же уже.
В среднем же длины связей оказываются сопоставимыми в диапазое 1.52-1.53 ангстрем, что соответствует действительности.
Построим распределение по потенциальной энергии¶
def plot_en_hist(short, long):
data_in = np.genfromtxt(short+'.xvg')
l = data_in[:,1]
fig, ax = plt.subplots(figsize=(10,5))
ax.hist(l, bins=250)
ax.set_title(f'Метод {long}')
ax.set_xlabel('Потенциальная энергия, кДж/моль')
for t, m in zip(types, methods):
plot_en_hist(t, m)
Формы практически всех распределений, кроме методов Берендсена и Андерсена, соответсвтуют распределению Больцамана с различными параметрами. Лучше всего ситуации при нормальных условиях соответствует метод Нуза-Хувера.
Время работы mdrun¶
Астрономическое время из логов в секундах:
Метод Берендсена: 4.132
Метод "Velocity rescale": 4.236
Метод Нуза-Хувера: 4.153
Метод Андерсена: 4.066
Метод стохастической молекулярной динамики: 4.666
Время вычисления каждым из методов оказалось сопоставимым. Однако быстрее всего сработал метод Андерсена, что не удивительно, поскольку он рассчитывает небольшие изменения системы, так как энергия меняется слабо. Для стохастической молекулярной динамики потребовалось больше всего времени.
В целом данные выглядят так, что при симуляции процессов, проходящих при комнатной температуре, лучше всего подходит метод Нуза-Хувера. Методы с большим разбросом по энергии будут оптимальны при расчете высоких температур.