pr7

Изучение работы методов контроля температуры в молекулярной динамике¶

Для начала построим файл топологии этана:

Скачали 5 конфигурационных файлов с разными алгоритмами контроля температуры. Получили из них tpr файлы (кроме an.mdp - с ним не работает). Потом для каждого из них запустили mdrun. Для анализа конвертировали результат (файл trr) в pdb. Посмотрим на то, что получилось, в PyMol.

Визуальный анализ¶

In [5]:

from IPython.display import Image

In [10]:

display('Метод Берендсена')
display(Image('et_be_movie.gif', format='png'))

'Метод Берендсена'

In [12]:

display('Метод "Velocity rescale"')
display(Image('et_vr_movie.gif', format='png'))

'Метод "Velocity rescale"'

In [13]:

display('Метод Нуза-Хувера')
display(Image('et_nh_movie.gif', format='png'))

'Метод Нуза-Хувера'

In [14]:

display('Метод стохастической молекулярной динамики')
display(Image('et_sd_movie.gif', format='png'))

'Метод стохастической молекулярной динамики'

Видим, что поведение этана разное (надеюсь, что все гифки отображаются нормально). Самым странным кажется метод стохастической молекулярной динамики - этан как-то странно сильно мечется.

Сравнение энергий¶

Теперь сравним потенциальную энегию и кинетическую энергию для всех систем.

In [16]:

import numpy as np
import matplotlib.pyplot as plt

In [25]:

systems = ['be', 'vr', 'nh', 'sd']
text = ['Берендсен', 'Velocity rescale', 'Нуз-Хувер', 'Стохастическая молекулярная динамика']

In [39]:

plt.subplots(2, 2, figsize=(20, 10))
for i, (n, txt) in enumerate(zip(data, text), 1):
    energy = np.loadtxt(f'et_{n}_en.xvg')
    t = energy[:,0]
    y_pot = energy[:,1]
    y_kin = energy[:,2]
    
    plt.subplot(2, 2, i)
    
    plt.plot(t, y_kin + y_pot, alpha=0.2)
    plt.plot(t, y_pot, alpha = 0.6)
    plt.plot(t, y_kin, alpha = 0.6)
    plt.title(txt)
    plt.legend(['Е полная', 'E потенциальная', 'E кинетическая'])
    plt.xlabel("Время", fontsize=10)
    plt.ylabel("Энергия",fontsize=10)
plt.show()

Видим, что только у Берендсена полная энергия поддерживается на одном уровне. В остальных системах - сильные колебания (самая нестабильная, кажется, Нуза-Хувера). А velocity rescale и стохастическая МД оказались очень похожи.

Распределение длины С-С связи¶

Сначала создали исходный индекс-файл с одной связью, и запустили утилиту distance.

In [55]:

plt.subplots(2, 2, figsize=(10, 8))
for i, (n, txt) in enumerate(zip(systems, text), 1):
    bond = np.loadtxt(f'bond_{n}.xvg')
    x = bond[:,0]
    y = bond[:,1]
    
    plt.subplot(2, 2, i)
    
    plt.title(txt)
    plt.bar(x, y, width=0.002)
    plt.xlim([0.11, 0.17])
plt.show()

Распределение потенциальной энергии¶

Распределение потенциальной энергии должно по идее иметь форму распределения Больцмана.

In [57]:

plt.subplots(2, 2, figsize=(10, 8))
for i, (n, txt) in enumerate(zip(systems, text), 1):
    energy = np.loadtxt(f'et_{n}_en.xvg')
    en_pot = energy[:,1]
    
    plt.subplot(2, 2, i)
    
    plt.hist(en_pot, bins=100)
    plt.title(txt)
plt.show()

Вроде, на распределение Больцмана больше похожи velocity rescaling и стохастическая МД.

Итак, по полученным результатам, можно сделать следующие выводы. Берендсен самый стабильный в плане энергии (и в плане распределения длин С-С связей тоже). Распределение Больцмана больше всех напоминают velocity rescale и стохастическая МД. Именно их я бы наверное и выбрала как самых реалистичных (они в целом похожи, но в стохастической МД есть движение молекулы как целого).

In [ ]: