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

Здесь мы на примере молекулы этана средствами программы GROMACS будем изучать различные методы контроля температуры в молекулярной динамике.

Подготовка файла координат и топологии

Файл координат берём с kodomo. Файл топологии допишем вручную на основе шаблона из задания. Файл доступен по ссылке. В файле учтены нужные типы атомов для силового поля OPLS и то, что тип функционала для торсионных углов в этом поле под номером 3.

Работа в GROMACS

Далее берём файлы конфигураций для молекулярной динамики, соответствующие методам контроля температуры Берендсена, Velocity rescale, Нуза — Хувера, стохастической молекулярной динамики. Для каждого метода применяем пайплайн из команд gmx grompp, gmx mdrun для получения небольшого фрагмента динамики.

Полученные результаты визуализируем в PyMol. Ниже приводятся анимации молекулярной динамики.

Берендсен

Velocity rescale

Нуз и Хувер

Стохастическая динамика

По виду можно сказать, что с балансом колебаний и вращений лучше всего справляется метод VR; в стохастической динамике вообще рвутся связи (??), что довольно странно. В термостате Берендсена малы колебания, у Нуза — Хувера нет вращения по торсионному углу.

Анализ распределений энергии и длин связей

Построим графики распределения потенциальной и кинетической энергии.

import numpy as np

be_energy = np.loadtxt("et_be_en.xvg")
be_energy

array([[ 0.0000000e+00, -7.9222590e+00,  8.7644000e-02],
       [ 1.0000000e-01, -4.3069930e+00,  5.6179330e+00],
       [ 2.0000000e-01,  7.4533410e+00,  1.0478281e+01],
       ...,
       [ 2.4980000e+02,  1.6541197e+01,  2.6601833e+01],
       [ 2.4990000e+02,  1.6722681e+01,  2.6362497e+01],
       [ 2.5000000e+02,  1.6646833e+01,  2.6528360e+01]])

t = be_energy[:, 0]
e_pot = be_energy[:, 1]
e_kin = be_energy[:, 2]

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(t, e_pot, lw=.5, label="$E_P$", alpha=.7)
plt.plot(t, e_kin, lw=.5, label="$E_{kin}$", alpha=.7)
plt.legend()
plt.title("Берендсен", fontsize=16)
plt.xlabel("Время", fontsize=14)
plt.ylabel("Энергия",fontsize=14)

Text(0, 0.5, 'Энергия')

vr_energy = np.loadtxt("et_vr_en.xvg")
t, e_pot, e_kin = (vr_energy[:, _] for _ in range(3))

plt.plot(t, e_pot, lw=.5, label="$E_P$", alpha=.7)
plt.plot(t, e_kin, lw=.5, label="$E_{kin}$", alpha=.7)
plt.legend()
plt.title("Velocity rescale", fontsize=16)
plt.xlabel("Время", fontsize=14)
plt.ylabel("Энергия",fontsize=14)

Text(0, 0.5, 'Энергия')

nh_energy = np.loadtxt("et_nh_en.xvg")
t, e_pot, e_kin = (nh_energy[:, _] for _ in range(3))

plt.plot(t, e_pot, lw=.5, label="$E_P$", alpha=.7)
plt.plot(t, e_kin, lw=.5, label="$E_{kin}$", alpha=.7)
plt.legend(loc="right")
plt.title("Нуз и Хувер", fontsize=16)
plt.xlabel("Время", fontsize=14)
plt.ylabel("Энергия",fontsize=14)

Text(0, 0.5, 'Энергия')

sd_energy = np.loadtxt("et_sd_en.xvg")
t, e_pot, e_kin = (sd_energy[:, _] for _ in range(3))

plt.plot(t, e_pot, lw=.5, label="$E_P$", alpha=.7)
plt.plot(t, e_kin, lw=.5, label="$E_{kin}$", alpha=.7)
plt.legend(loc="right")
plt.title("Стохастическая динамика", fontsize=16)
plt.xlabel("Время", fontsize=14)
plt.ylabel("Энергия",fontsize=14)

Text(0, 0.5, 'Энергия')

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

paths = ["et_be_en.xvg",
         "et_vr_en.xvg",
         "et_nh_en.xvg", 
         "et_sd_en.xvg",]
titles = ["Берендсен",
          "Velocity rescale",
          "Нуз и Хувер",
          "Стохастическая динамика"]

def plot(which):
    energy = np.loadtxt(paths[which])
    t, e_pot, e_kin = (energy[:, _] for _ in range(3))
    e_pot -= e_pot[0]
    
    plt.plot(t, e_pot, lw=.5, label="$E_P$", alpha=.7)
    plt.plot(t, e_kin, lw=.5, label="$E_{kin}$", alpha=.7)
    plt.legend(loc="right")
    plt.title(titles[which], fontsize=16)
    plt.xlabel("Время", fontsize=14)
    plt.ylabel("Энергия",fontsize=14)

plot(0)

plot(1)

plot(2)

plot(3)

В термостате Берендсена виды энергии практически не переходят друг в друга (отклонения малы). В случае Нуза — Хувера у энергии очень большой разброс, но в целом она тяготеет к околонулевым значениям. Графики VR и стохастической динамики довольно схожи, выглядят неплохо, там характерное отклонение имеет тот же порядок, что и величина энергии, поэтому виды энергии эффективно друг в друга переходят.

Займёмся длинами связей. Команда из задания сразу задаёт гистограмму распределения. Надо только правильно прочитать и визуализировать файл.

bond_be = np.loadtxt("bond_be.xvg")
existing_bonds_be = bond_be[bond_be[:,1]>0]
plt.bar(existing_bonds_be[:,0], existing_bonds_be[:,1],width=.001)
plt.title("Гистограмма распределения длин связей \nдля термостата Берендсена", fontsize=16)
plt.xlabel("Длина связи", fontsize=14)
plt.ylabel("Плотность", fontsize=14)

Text(0, 0.5, 'Плотность')

bond_vr = np.loadtxt("bond_vr.xvg")
existing_bonds_vr = bond_vr[bond_vr[:,1]>0]
plt.bar(existing_bonds_vr[:,0], existing_bonds_vr[:,1],width=.001)
plt.title("Гистограмма распределения длин связей \nдля термостата Velocity rescale", fontsize=16)
plt.xlabel("Длина связи", fontsize=14)
plt.ylabel("Плотность", fontsize=14)

Text(0, 0.5, 'Плотность')

bond_nh = np.loadtxt("bond_nh.xvg")
existing_bonds_nh = bond_nh[bond_nh[:,1]>0]
plt.bar(existing_bonds_nh[:,0], existing_bonds_nh[:,1],width=.001)
plt.title("Гистограмма распределения длин связей \nдля термостата Нуза — Хувера", fontsize=16)
plt.xlabel("Длина связи", fontsize=14)
plt.ylabel("Плотность", fontsize=14)

Text(0, 0.5, 'Плотность')

bond_sd = np.loadtxt("bond_sd.xvg")
existing_bonds_sd = bond_sd[bond_sd[:,1]>0]
plt.bar(existing_bonds_sd[:,0], existing_bonds_sd[:,1],width=.001)
plt.title("Гистограмма распределения длин связей \nдля термостата стохастической динамики", fontsize=16)
plt.xlabel("Длина связи", fontsize=14)
plt.ylabel("Плотность", fontsize=14)

Text(0, 0.5, 'Плотность')

Здесь видно, что для термостата Нуза — Хувера распределение крайне асимметрично, поэтому он работает плохо. Для термостата Берендсена наблюдается очень маленькая дисперсия (заморозка колебаний). Наконец, в случае VR и SD всё хорошо. Теперь стало понятно, что никакие связи в стохастической динамике не рвались, это просто молекула выходила за границы ячейки и транслировалась с другой стороны.

Проверим также соответствие распределения потенциальной энергии по времени для частицы распределению Больцмана.

def plot_distribution(which):
    energy = np.loadtxt(paths[which])
    t, e_pot, e_kin = (energy[:, _] for _ in range(3))
    
    plt.hist(e_pot, bins=20)
    plt.title(titles[which], fontsize=16)
    
    if which > 1:
        plt.xlabel("Потенциальная энергия", fontsize=14)
    plt.ylabel("Плотность вероятности",fontsize=14)

plt.subplots(2, 2, figsize=(10,8))

for i in range(4):
    plt.subplot(2, 2, i + 1)
    plot_distribution(i)

Распределение Больцмана несимметричное, с тяжёлым правым хвостом. Лучше всего ему соответствуют VR и SD, хотя тоже неидеально (мода немного не совпадает с минимумом значения).

Быстродействие

Ниже представлена таблица времени работы (wall time) всех алгоритмов.

Алгоритм	Время, с
Be	22.0
VR	22.2
NH	23.1
SD	21.4

Заключение

Исследовав 4 термостата, можем заключить, что лучшим образом себя показали VR и SD. К тому же термостат SD отработал быстрее всего, что также говорит в его пользу.