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

Подготовим все файлы

Напишем файл топологии молекулы

et=open('et.top','w')
et.write('''#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_135   1    ETH      C1      1    -0.189      12.01
    2   opls_135   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
    1     2     6     1
    1     2     7     1
    1     2     8     1
    6     2     7     1
    6     2     8     1
    7     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 ]
; list of atoms 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''')
et.close()

Подготовим команды для GROMACS, которые затем будем использовать в стороннем ресурсе.

При этом сразу "создадим" команду для рассчета длины связи.

При этом учтем, что у нас имеются 5 конфигурационных файлов с разными алгоритмами контроля температуры:

• be.mdp - метод Берендсена для контроля температуры.
• vr.mdp - метод "Velocity rescale" для контроля температуры.
• nh.mdp - метод Нуза-Хувера для контроля температуры.
• an.mdp - метод Андерсена для контроля температуры.
• sd.mdp - метод стохастической молекулярной динамики.

types = ['be', 'vr', 'nh', 'an', 'sd']

for i in types:
    print(f'gmx grompp -f {i}.mdp -c etane.gro -p et.top -o et_{i}.tpr')
    print(f'gmx mdrun -deffnm et_{i} -v -nt 1')
    print(f'gmx trjconv -f et_{i}.trr -s et_{i}.tpr -o et_{i}.pdb')
    print(f'gmx energy -f et_{i}.edr -o et_{i}_en.xvg -xvg none')
    print(f'gmx distance -f et_{i}.trr -s et_{i}.tpr -oh bond_{i}.xvg -n b.ndx -xvg none -select 0')
    print('--------------------------------------------')

gmx grompp -f be.mdp -c etane.gro -p et.top -o et_be.tpr
gmx mdrun -deffnm et_be -v -nt 1
gmx trjconv -f et_be.trr -s et_be.tpr -o et_be.pdb
gmx energy -f et_be.edr -o et_be_en.xvg -xvg none
gmx distance -f et_be.trr -s et_be.tpr -oh bond_be.xvg -n b.ndx -xvg none -select 0
--------------------------------------------
gmx grompp -f vr.mdp -c etane.gro -p et.top -o et_vr.tpr
gmx mdrun -deffnm et_vr -v -nt 1
gmx trjconv -f et_vr.trr -s et_vr.tpr -o et_vr.pdb
gmx energy -f et_vr.edr -o et_vr_en.xvg -xvg none
gmx distance -f et_vr.trr -s et_vr.tpr -oh bond_vr.xvg -n b.ndx -xvg none -select 0
--------------------------------------------
gmx grompp -f nh.mdp -c etane.gro -p et.top -o et_nh.tpr
gmx mdrun -deffnm et_nh -v -nt 1
gmx trjconv -f et_nh.trr -s et_nh.tpr -o et_nh.pdb
gmx energy -f et_nh.edr -o et_nh_en.xvg -xvg none
gmx distance -f et_nh.trr -s et_nh.tpr -oh bond_nh.xvg -n b.ndx -xvg none -select 0
--------------------------------------------
gmx grompp -f an.mdp -c etane.gro -p et.top -o et_an.tpr
gmx mdrun -deffnm et_an -v -nt 1
gmx trjconv -f et_an.trr -s et_an.tpr -o et_an.pdb
gmx energy -f et_an.edr -o et_an_en.xvg -xvg none
gmx distance -f et_an.trr -s et_an.tpr -oh bond_an.xvg -n b.ndx -xvg none -select 0
--------------------------------------------
gmx grompp -f sd.mdp -c etane.gro -p et.top -o et_sd.tpr
gmx mdrun -deffnm et_sd -v -nt 1
gmx trjconv -f et_sd.trr -s et_sd.tpr -o et_sd.pdb
gmx energy -f et_sd.edr -o et_sd_en.xvg -xvg none
gmx distance -f et_sd.trr -s et_sd.tpr -oh bond_sd.xvg -n b.ndx -xvg none -select 0
--------------------------------------------

Дополнительные комментарии:

• для метода Берендсена пришлось указатьмаксимальное количество предупреждений, так как программа предупреждала, что данный термостат не создает правильного распределения кинетической энергии.
• для метода Андерсена было необходимо изменить файл с настройками, так как для этого метода необходим интегратор md-vv, а не просто md

Визуализация в PyMol

import imageio
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import os

import pymol
pymol.finish_launching()
from pymol import cmd,stored
from IPython.display import Image, display

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

cmd.do(f'''
reinitialize everything
load et_be.pdb
mset 1-{cmd.count_states()}
set ray_trace_frames,1
set cache_frames,0
bg_color white
set ray_opaque_background, on
mclear
mpng ./et_be/mov
''')

filenames = sorted(os.listdir('./et_be/'))[1:]
with imageio.get_writer('et_be.gif', mode='I', fps=20,loop=0) as writer:   
    for filename in filenames:
        filepath = os.path.join('./et_be', filename)
        image = imageio.imread(filepath)
        writer.append_data(image)

display(Image(filename='et_be.gif', format='png'))

Метод "Velocity rescale"

cmd.do(f'''
reinitialize everything
load et_vr.pdb
mset 1-{cmd.count_states()}
set ray_trace_frames,1
set cache_frames,0
bg_color white
set ray_opaque_background, on
mclear
mpng ./et_vr/mov
''')

filenames = sorted(os.listdir('./et_vr/'))[1:]
with imageio.get_writer('et_vr.gif', mode='I', fps=20,loop=0) as writer:   
    for filename in filenames:
        filepath = os.path.join('./et_vr', filename)
        image = imageio.imread(filepath)
        writer.append_data(image)

display(Image(filename='et_vr.gif', format='png'))

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

cmd.do(f'''
reinitialize everything
load et_nh.pdb
mset 1-{cmd.count_states()}
set ray_trace_frames,1
set cache_frames,0
bg_color white
set ray_opaque_background, on
mclear
mpng ./et_nh/mov
''')

filenames = sorted(os.listdir('./et_nh/'))[1:]
with imageio.get_writer('et_nh.gif', mode='I', fps=20,loop=0) as writer:   
    for filename in filenames:
        filepath = os.path.join('./et_nh', filename)
        image = imageio.imread(filepath)
        writer.append_data(image)

display(Image(filename='et_nh.gif', format='png'))

Метод Андерсена

cmd.do(f'''
reinitialize everything
load et_an.pdb
mset 1-{cmd.count_states()}
set ray_trace_frames,1
set cache_frames,0
bg_color white
set ray_opaque_background, on
mclear
mpng ./et_an/mov
''')

filenames = sorted(os.listdir('./et_an/'))[1:]
with imageio.get_writer('et_an.gif', mode='I', fps=20,loop=0) as writer:   
    for filename in filenames:
        filepath = os.path.join('./et_an', filename)
        image = imageio.imread(filepath)
        writer.append_data(image)

display(Image(filename='et_an.gif', format='png'))

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

cmd.do(f'''
reinitialize everything
load et_sd.pdb
mset 1-{cmd.count_states()}
set ray_trace_frames,1
set cache_frames,0
bg_color white
set ray_opaque_background, on
mclear
mpng ./et_sd/mov
''')

filenames = sorted(os.listdir('./et_sd/'))[1:]
with imageio.get_writer('et_sd.gif', mode='I', fps=20,loop=0) as writer:   
    for filename in filenames:
        filepath = os.path.join('./et_sd', filename)
        image = imageio.imread(filepath)
        writer.append_data(image)

display(Image(filename='et_sd.gif', format='png'))

Вывод

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

Сравнение энергий для каждой из систем

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

import numpy as np
import matplotlib.pyplot as plt
def plot_energy(short, long):
    data_in = np.genfromtxt(f'et_{short}_en.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)

Интересно добавить на графики еще и Полную энергию системы, чтобы посмотреть еще и на то, как связаны Кинетическая, Потенциальная и Полная энергии с течением времени в каждой из систем.

import numpy as np
import matplotlib.pyplot as plt
def plot_energy(short, long):
    data_in = np.genfromtxt(f'et_{short}_en.xvg')
    x, pot, kin, full = data_in[:,0], data_in[:,1], data_in[:,2], data_in[:,3]
    
    fig, ax = plt.subplots(figsize=(10,5))
    ax.plot(x, pot, label='Потенциальная E')
    ax.plot(x, kin, label='Кинетическая E')
    ax.plot(x, full, 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)

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

Также заметим то, что для стохастической молекулярной динамики и метода Андерсена графики отражают чуть ли не одно и то же.

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

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

def plot_bonds(short, long):
    data_in = np.genfromtxt(f'bond_{short}.xvg')
    x, l = data_in[:,0], data_in[:,1]
    
    fig, ax = plt.subplots(figsize=(10,5))
    ax.bar(x, l, width=0.001)
    ax.set_title(f'Метод {long}')
    ax.set_xlabel('Длина связи, нм')

methods = ['Берендсена', '"Velocity rescale"', 'Нуза-Хувера',
           'Андерсена', 'стохастической молекулярной динамики']

for t, m in zip(types, methods):
    plot_bonds(t, m)

В среднем же длины связей оказываются сопоставимыми в диапазое 1.52-1.53 ангстрем, что соответствует действительности.

В целом каждый из методов имеет четкий пик в указанном диапазоне. Что-то странное видно только в методе Нуза-Хувера: есть какой-то странный кластер из связей длиной 1.13-1.17 ангстрем.

Время работы mdrun и общий вывод

Астрономическое время из логов в секундах:

• Метод Берендсена: 9.759
• Метод "Velocity rescale": 9.999
• Метод Нуза-Хувера: 9.641
• Метод Андерсена: 10.125
• Метод стохастической молекулярной динамики: 9.689

В целом, время работы для каждого из методов оказалось вполне сопоставимым. Поэтому особо нечего сказать про различия во времени работы.

Исходя из всех данных, представленных в этой работе, можно сказать следующее:

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