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

Эта работа была посвещена изучению того, как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования это одна молекула этана. Для нее вручную был дополнен файл топологии et.top.

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

Необходимо было запустить молекулярную динамику GROMACS для этана с разными параметрами контроля температуры и по результатам сравнить данные 5 методов. Для этого было сделано следующее:
1) Построены входные файлы для молекулярно-динамического движка mdrun с помощью grompp;
2) Получили 5 tpr файлов. Теперь для каждого из них запустили mdrun для оптимизации геометрии молекулы;
3) Для каждой из 5 систем провели конвертацию в pdb для визульного анализа.

In [13]:
import subprocess, os
my_env = os.environ.copy()
for i in ["be","vr","nh","an","sd"]:
    do_grompp="grompp -f %s.mdp -c et.gro -p et.top -o et_%s.tpr" %(i,i)
    do_mdrun="mdrun -deffnm et_%s -v -nt 1" %i
    make_pdb="echo 1 | trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb" %(i,i,i)
    p1 = subprocess.Popen(do_grompp, 
                     shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env) 
    p2 = subprocess.Popen(do_mdrun, 
                     shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)
    p3 = subprocess.Popen(make_pdb, 
                     shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)

В результате было получено 5 pdb файлов: et_be.pdb, et_vr.pdb, et_nh.pdb, et_an.pdb, et_sd.pdb. Все эти файлы были проанализированы в pymol. В результате можно сделать следующие выводы:

  • При контроле температуры методом Берендсена молекула этана сначала быстро вращается вокруг СС-связи, начинает вращаеться вокруг других осей, совершая разворот молекулы на 360° и немного меняя уклон относительно горизонта, потом вращательное движение вокруг С-С связи прекращается и переходит в колебальное движение водородов относительно С-С связи;
  • При контроле температуры методом "Velocity rescale" молекула этана совершает вращатеные или колебательные движения вокруг С-С связи и немного разворачивается (примерно на 15°) в горизонтальной плоскости;
  • При контроле температуры методом Нуза-Хувера молекула этана совершает вращательные или колебательные движения вокруг СС-связи, при этом молекула целиком остается на месте;
  • При контроле температуры методом Андерсена молекула этана совершает только небольшие колебательные движения, в ходе которых немного меняется длинна С-С связи;
  • При контроле температуры методом стохастической молекулярной динамики молекула этана быстро и хаотично вращается вокруг различных "случайных" осей.

4) Далее рассчитали потенциальную и кинетическую энергии для каждой из 5 систем.

In [15]:
import subprocess, os
my_env = os.environ.copy()
for i in ["be","vr","nh","an","sd"]:
    pot_energy="echo 10 | g_energy -f et_%s.edr -o et_%s_en_pot.xvg -xvg none" %(i,i)
    kin_energy="echo 11 | g_energy -f et_%s.edr -o et_%s_en_kin.xvg -xvg none" %(i,i)
    p4 = subprocess.Popen(pot_energy, 
                     shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)
    p5 = subprocess.Popen(kin_energy, 
                     shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)

Получили 5 файлов с информацией о потенциальной энергии молекулы этана (et_vr_en_pot.xvg, et_an_en_pot.xvg, et_be_en_pot.xvg, et_nh_en_pot.xvg, et_sd_en_pot.xvg) и 5 файлов с информацией о кинетической энергии (et_vr_en_kin.xvg, et_an_en_kin.xvg, et_be_en_kin.xvg, et_nh_en_kin.xvg, et_sd_en_kin.xvg).

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

In [20]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
titles = {
'an': 'Anderson method of temperature control',
'be': 'Berendsen method of temperature control',
'nh': 'Nose-Hoover method of temperature control',
'vr': '"Velocity rescale" method of temperature control',
'sd': 'Stochastic molecular dynamics method'}

for i in ["be","vr","nh","an","sd"]:
    a_pot = np.loadtxt('et_%s_en_pot.xvg' % i)
    t_pot = a_pot[:,0]
    y_pot = a_pot[:,1]
    a_kin = np.loadtxt('et_%s_en_kin.xvg' % i)
    t_kin = a_kin[:,0]
    y_kin = a_kin[:,1]
    plt.plot(t_pot, y_pot, "bo", t_kin, y_kin, "ro", alpha=0.6)
    red_patch = mpatches.Patch(color='red', label='Kinetic energy')
    blue_patch = mpatches.Patch(color='blue', label='Potential energy')
    plt.legend(handles=[red_patch,blue_patch])
    plt.title(titles[i]) 
    plt.xlabel("Time (ps)")
    plt.ylabel("Energy (kJ/mol)")
    plt.show()

Во всех случаях кинетическая и понетциальные энергии практически совпадают. Метод "Velocity rescale" и метод стохастической молекулярной динамики очень похожи друг на друга, диапозон значений внутренней энергии системы от 0 до 60 кДж/моль. В случае контроля темтературы методом Андерсена внутренняя энергия системы близка к нулю (диапозон от 0 до 1,4 кДж/моль). При контроле температуры методом Нуза-Хувера внутренняя энергия системы в целом выше, чем при контроле температуры другими методами, а также наблюдаются аномально высокие значения кинетической энергии (>500 кДж/моль. При контроле температуры методом Берендсена потенциальная и кинетическая энергии находятся в узком диапазоне, который со временем уменьшается: в первый момент времени диапозон 0-40 кДж/моль, потом 20-35 кДж/моль и ближе к концу 23-27 кДж/моль.

6)Далее рассмотрели длинны связи С-С за время моделирования.

In [16]:
for i in ["be","vr","nh","an","sd"]:
    do_bond = 'g_bond -f et_%s.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none' %(i,i,i)
    p6 = subprocess.Popen(do_bond, 
                     shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)

Получили 5 файлов с информацией о рассчитаной длине С-С связи для каждого метода: bond_an.xvg, bond_be.xvg, bond_vr.xvg, bond_nh.xvg, bond_sd.xvg.

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

In [22]:
for i in ["be","vr","nh","an","sd"]:
    a = np.loadtxt("bond_%s.xvg" % i)
    x = a[:,0]
    y = a[:,1]
    width = 0.0003
    plt.bar(x, y, width, color="blue")
    plt.title(titles[i])
    plt.xlabel("CC-bond length (nm)")
    plt.show()

Распеределения длин С-С связей наиболее похожи на нормальные при использовании метода Берендсена, метода Нуза-Хувера и метода стохастической молекулярной динамики, максимумы данных распеределение попадают на 1.53Å, что соответстует дейсвительности. Другие два графика не похожи на нормальное распределение и имеют максимумы в других местах.

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

  • метод Андерсена изменяет только длины С-С связи, имеет значения энергий молекулы близкие к нулю, поэтому этот метод не подходит для моделирования системы из одной молекулы;
  • метод Берендсена делает хорошее распределение длин связей, но внутренняя энергия системы находится в маленьком диапозоне и этот диапозон уменьшается со временем, что выглядит мало вероятным для реальной системы;
  • метод Нуза-Хувера дает сильные выбросы в значениях кинетической энергии молекулы, поэтому этот метод тоже не подходит для моделирования системы из одной молекулы;
  • метод "Velocity rescale" дает не нормальное распределение длин связей с максимум не в 1,53Å, и этот метод тоже не подходит;
  • метод стохастической молекулярной динамики имеет хорошее распределение связей, хорошие значения кинетической и потенциальной энергий в молекуле, но при визуализации динамики молекулы мы видим хаотичные изменения положения молекулы, но все же я считаю, что этот метод самый лучший для реализации контроля температуры в молекулярной динамике одной молекулы.