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

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

Изменяем файл с топологией:

  • типы атомов нужно изменить: CX на opls_139, HX - на opls_140
  • дописываем связи и углы
  • тип функционала для торсионных углов меняем на 3
In [3]:
#смотрим файл с топологией
for line in open('et.top', encoding='utf-8').readlines():
    print(line.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_139      1    ETH      C1      1    -0.189      12.01
2   opls_139      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
6     2     7     1
7     2     8     1
6     2     8     1
1     2     6     1
1     2     7     1
1     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 ]
; список атомов 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

Скачиваем все необходимые файлы. Далее, с помощью bash скриптов я:

  • построила входные файлы .tpr для mdrun с помощью grompp
  • для каждого из файлов запустила mdrun
  • конвертировала файлы с результатами в pdb (trjconv)
  • сравнила потенциальную энергию и кинетическую (g_energy)
  • провела анализ длин связей при помощи g_bond

Сначала визуально оценим то, что получилось, в Pymol

  • Андерсен: молекула почти не двигается, длины связей и углы немножко колеблются вокруг одного положения. Сохраняется идеальная заторможенная конформация.
  • Берендсен: происходит постоянное вращение вокруг C-C связи. Молекула при этом всё время сохраняет заторможенную конформацию, то есть смена одной заторможенной конформации на другую мгновенная. Также, молекула, видимо, приобретает дополнительное ускорение и целиком подворачивается в одном направлении.
  • Velocity-rescale: вращение более спокойное, не всегда в одну сторону. Молекула бывает в заслонённой конформации, но чаще в заторможенной.
  • Нуз-Хувер: движение похоже на Velocity-rescale, сама молекула никуда не смещается.
  • Стохастическая МД: вся молекула хаотично поворачивается,при этом остаётся в заторможенной конформации.

Исходя из визуального анализа, ближе всего к правде Velocity-rescale и Нуз-Хувер.

Построим графики изменения энергий

In [4]:
import numpy as np
import matplotlib.pyplot as plt
In [16]:
algs = ['be', 'vr', 'nh', 'an', 'sd']
names = {'be': 'Berendsen', 'vr': 'Velocity-rescale', 'nh':'Nooze-Hoover', 'an':'Anderson', 'sd':'Stochastic MD'}
In [17]:
fig, axs = plt.subplots(5, figsize=(10, 18))
for i, d in enumerate(algs):
    energy = np.loadtxt(f'pr7/et_{d}.en.xvg')
    t = energy[:,0]
    e_pot = energy[:,1]
    e_kin = energy[:,2]
    axs[i].plot(t, e_pot, c='blue', alpha=0.6, label='Potential energy')
    axs[i].plot(t, e_kin, c='red', alpha=0.6, label='Kinetic energy')
    axs[i].legend()
    axs[i].set_title(names[algs[i]])
plt.show()

В случае Берендсена обе энергии быстро стабилизируются. Нуз-Хувер - резкие скачки кинетической энергии, где-то в середине энергии тоже стабилизируются. Остальные три метода: постоянные флуктуации обеих энергий.

Построим графики распределения длин связей.

In [24]:
fig, axs = plt.subplots(5, figsize=(10, 18))
for i, d in enumerate(algs):
    x = np.loadtxt(f'pr7/bond_{d}.xvg')
    axs[i].bar(x[:,0], x[:,1], width=0.0005)
    axs[i].set_title(names[algs[i]])
    axs[i].set_xlim(0.14, 0.163)
plt.show()

Дисперсия распределений больше всего у стохастической МД, у Velocity-rescale тоже большая. У Андерсона высокая плотность посередине и по краям: не очень реалистично. Распределения для Берендсена и Нуза-Хувера похожи, у последнего побольше дисперсия.

Построим распределения по энергиям.

In [25]:
fig, axs = plt.subplots(5, figsize=(10, 18))
for i, d in enumerate(algs):
    energy = np.loadtxt(f'pr7/et_{d}.en.xvg')
    t = energy[:,0]
    e_pot = energy[:,1]
    e_kin = energy[:,2]
    axs[i].hist(e_pot + e_kin, bins=100, color='purple')
    axs[i].set_title(names[algs[i]])
plt.show()

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

P.S. Скрипты и все файлы можно найти в /home/students/y17/darya_by/MM/pr7 на kodomo

In [ ]: