Изучение работы методов контроля температуры в молекулярной динамике¶
В данном практикуме рассматривается контроль температуры в рамках программы GROMACS. В процессе симуляции система рассматривается как микроканонический ансамбль, c постоянной внутренней энергией. Это означет, что нужен дополнительный контроль температуры, чтобы система не разогревалась и не замерзала из-за перехода потенциальной энергии в кинетическую и наоборот.
Подготовка файлов координат и топологии¶
Я взяла готовые файлы координат и топологии для силового поля OPLS.
Для торсионных углов я указала тип потенциала "3" -- потенциал Рюкарта-Бельманса, который представляет собой полином из косинусов.
Типы атомов для углерода и водорода соответственно -- opls_135 и opls_140.
Также я дописала недостающие строки для связей и углов.
Запуск динамики¶
Для генерации .trp файлов, запуска mdrun, генерации pdb файлов, вывода кинетической/потенциальной/полной энергий, получения длин CC связей я написала bash-скрипт, который проходится по mdp файлам для разных термостатов:
#!/bin/bash
shopt -s nullglob
for file in *.mdp; do
basename="${file%.mdp}"
gmx grompp -f "$file" -c et.gro -p et.top -o "et_${basename}.tpr"
gmx mdrun -deffnm "et_${basename}" -v -nt 2
echo 0 | gmx trjconv -f "et_${basename}.trr" -s "et_${basename}.tpr" -o "et_${basename}.pdb"
echo "8 9 10 0" | gmx energy -f "et_${basename}.edr" -o "et_${basename}_en.xvg" -xvg none
echo 0 | gmx distance -f "et_${basename}.trr" -s "et_${basename}.tpr" -oh "bond_${basename}.xvg" -n b.ndx -xvg none
done
echo "done"
Визуальный анализ¶
Метод Берендсена¶
На вид молекула просто крутиться как руль корабля:
from IPython.display import Image
Image(url='be.gif')
Метод Velocity rescale¶
Image(url='vr.gif')
Метод Нуза-Хувера¶
Торсионные углы особо не меняются :(
Image(url='nh.gif')
Метод стохастической молекулярной динамики¶
Этан оказался слишком подвижным и убегающим. Но это не является недостатком.
Image(url='sd.gif')
На вид наиболее привлекательным методом кажется velocity rescale. И торсионные углы меняются, и длина CC связи. Результат стохастической молекулярной динамики также кажется довольно убедительным.
Анализ энергий¶
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
s = {"be":"метод Берендсена",
"vr": "Velocity rescale",
"nh": "метод Нуза-Хувера",
"sd": "метод стохастической молекулярной динамики",
"an": "метод Андерсена",
}
fig, axes = plt.subplots(nrows=4, ncols=1,
figsize=(10, 2.5 * 4),
sharex=True,
dpi=100)
for (suff, name), ax in zip(s.items(), axes):
pth = Path(f"./et_{suff}_en.xvg")
if pth.exists():
data = np.loadtxt(pth, comments=['@', '#'])
time = data[:, 0]
ax.plot(time,
data[:, 1],
label=f'{name} потенциальная энергия',
linewidth=1)
ax.plot(time,
data[:, 2],
label=f'{name} кинетическая энергия',
linewidth=1)
ax.grid(True, linestyle='--', alpha=0.5)
ax.legend(loc='upper right', fontsize=10)
ax.set_ylabel('Энергия', fontsize=12)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
axes[-1].set_xlabel('Время (ps)', fontsize=14)
axes[0].set_title('Сравнение энергий для различных термостатов', fontsize=16, fontweight='bold', pad=15)
plt.tight_layout()
И еще посмотрим на гистограммы распределения энергий:
s = {"be":"метод Берендсена",
"vr": "Velocity rescale",
"nh": "метод Нуза-Хувера",
"sd": "метод стохастической молекулярной динамики",
"an": "метод Андерсена",
}
fig, axes = plt.subplots(nrows=4, ncols=1,
figsize=(10, 2.5 * 4),
sharex=True,
dpi=100)
for (suff, name), ax in zip(s.items(), axes):
pth = Path(f"./et_{suff}_en.xvg")
if pth.exists():
data = np.loadtxt(pth, comments=['@', '#'])
time = data[:, 0]
ax.hist(
data[:, 1],
label=f'{name} потенциальная энергия',
linewidth=1)
ax.hist(
data[:, 2],
label=f'{name} кинетическая энергия',
linewidth=1)
ax.grid(True, linestyle='--', alpha=0.5)
ax.legend(loc='upper right', fontsize=10)
ax.set_ylabel('', fontsize=12)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
axes[-1].set_xlabel('Энергия', fontsize=14)
axes[0].set_title('Распределение энергий для различных термостатов', fontsize=16, fontweight='bold', pad=15)
plt.tight_layout()
Уже всего интервал колебаний энергий получился для метода Берендсена, однако хорошо ли это? Колокол сильно заужен, что означает, что система скорее всего не исследует все микросостояния. Методы velocity rescale и стохастической молекулярной динамики выглядят красиво и приятно, похоже на гауссиану, что соотвествует распределению Больцмана. Метод Нуза-Хувера выглядит не очень похоже на гауссиану.
Анализ расстояний CC¶
s = {"be":"метод Берендсена",
"vr": "Velocity rescale",
"nh": "метод Нуза-Хувера",
"sd": "метод стохастической молекулярной динамики",
"an": "метод Андерсена",
}
fig, axes = plt.subplots(nrows=4, ncols=1,
figsize=(10, 2.5 * 4),
sharex=True,
dpi=100)
for (suff, name), ax in zip(s.items(), axes):
pth = Path(f"./bond_{suff}.xvg")
if pth.exists():
data = np.loadtxt(pth, comments=['@', '#'])
time = data[:, 0]
ax.plot(time,
data[:, 1],
label=f'{name}',
linewidth=1)
ax.grid(True, linestyle='--', alpha=0.5)
ax.legend(loc='upper right', fontsize=10)
ax.set_ylabel('', fontsize=12)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
axes[-1].set_xlabel('Расстояние CC', fontsize=14)
axes[0].set_title('Гистограммы распределений длин CC', fontsize=16, fontweight='bold', pad=15)
plt.tight_layout()
Уже всего пик получился для метода Берендсена (в районе 0.15 А), алгоритм искусственно обрезает флуктуации. Более разлапистые гауссианы получились для остальных методов.
Быстродействие термостатов¶
!grep -H "Performance:" *.log
et_be.log:Performance: 7059.920 0.003 et_mdout.log:Performance: 6988.325 0.003 et_nh.log:Performance: 3333.883 0.007 et_sd.log:Performance: 6970.993 0.003 et_vr.log:Performance: 7067.040 0.003
Во второй колонке находсятся данные о быстродействие алгоритмов: сколько наносекунд симуляции в день будет считаться тот или иной метод. Из этих данных можно сделать вывод, что наиболее быстро справляются методы Берендсена и velocity rescale, самый медленный метод -- Нуза-Хувера.
mdout случайно добавился при итерированию по mdp файлам папки
Выводы¶
Лучше всего соотвествуют распределению Больцмана методы velocity rescale и метод стозастической динамики. Также эти методы визуально наиболее реалистичны. Самыми быстрыми методами являются методы Берендмена и velocity rescale. Таким образом наилучшим методом из рассмаотренных вероятно является velocity rescale.
метод Андерсена содержит какую-то ошибку в mdp файле, поэтому я его пропустила