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

Данный практикум был посвящен тому, как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования - молекула этана.

Нам был дан файл координат этана (etane.gro) и файл его топологии (et.top). Файл топологии был исправлен вручную (были добавлены недостающие связи и углы и добавлены типы атомов из файла /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp).

Также нам были даны 5 файлов с разными параметрами контроля температуры:

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

Необходимо было запустить молекулярную динамику GROMACS для этана с разными параметрами контроля температуры и по результатам сравнить данные 5 методов.

In [2]:
#запускаем динамику с разными параметрами для контроля температуры и получаем соответствующие pdb-файлы
import subprocess
my_env = os.environ.copy()
for method in ["be","vr","nh","an","sd"]:
    #Строим входные файлы для молекулярно-динамического движка mdrun с помощью grompp
    command1 = 'grompp -f %s.mdp -c etane.gro -p et.top -o et_%s.tpr' %(method,method)
    p1 = subprocess.Popen(command1,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)
In [4]:
import time
for method in ["be","vr","nh","an","sd"]:
    #Запускаем mdrun и сразу считаем время работы для каждого алгоритма
    command2 = 'mdrun -deffnm et_%s -v -nt 1' %method
    start_time = time.time()
    subprocess.call(command2,shell=True)
    print("Execution time for %s-method: %s seconds" % (method,time.time() - start_time))

Execution time for be-method: 4.28529787064 seconds
Execution time for vr-method: 4.29768586159 seconds
Execution time for nh-method: 4.37534213066 seconds
Execution time for an-method: 4.21546292305 seconds
Execution time for sd-method: 4.88758587837 seconds

Сразу можно отметить, что время быстродействия алгоритмов отличается не сильно. Самыми быстрыми можно считать алгоритмы an, be и vr, алгоритм nh немного помедленнее и алгоритм sd медленнее всех.

In [14]:
for method in ["be","vr","nh","an","sd"]:
    #Проводим конвертацию выходных файлов в pdb
    command3 = 'echo 1 | trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb' %(method,method,method)
    p3 = subprocess.Popen(command3,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)

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

Для каждого метода скачали полученный pdb и посмотрели в PyMol (на анимациях ниже показаны только 35 первых кадров)

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

Видно, что колебания молекулы минимальные, меняется только длина С-С связи

In [17]:
from IPython.display import Image
Image(url='an_movie.gif')
Out[17]:

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

Довольно быстрые вращения вокруг С-С связи вместе с небольшими изменениями длины С-С связи

In [18]:
Image(url='be_movie.gif')
Out[18]:

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

Немного похож на предыдущий, но амплитуда и скорость поворотов вокруг С-С связи меньше

In [19]:
Image(url='nh_movie.gif')
Out[19]:

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

Быстрые хаотичные движения молекулы целиком

In [20]:
Image(url='sd_movie.gif')
Out[20]:

метод "Velocity rescale"

По характеру движения похож на nh и be, небольшие и менее резкие повороты вокруг С-С связи

In [21]:
Image(url='vr_movie.gif')
Out[21]:

При визуальном анализе движения молекулы этана видно, что наиболее правдоподобными кажутся методы nh, vr и be, наименее правдоподобными - an и особенно sd.

Сравнение потенциальной и кинетической энергии систем

In [31]:
for method in ["be","vr","nh","an","sd"]:
    #Потенциальная энергия
    command1 = 'echo 10 | g_energy -f et_%s.edr -o et_%s_pot_en.xvg -xvg none' %(method,method)
    p1 = subprocess.Popen(command1,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)
    #Кинетическая энергия
    command2 = 'echo 11 | g_energy -f et_%s.edr -o et_%s_kin_en.xvg -xvg none' %(method,method)
    p2 = subprocess.Popen(command2,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=my_env)
In [1]:
#Строим графики
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
for method in ["be","vr","nh","an","sd"]:
    a_pot = np.loadtxt('et_%s_pot_en.xvg' % method)
    t_pot = a_pot[:,0]
    y_pot = a_pot[:,1]
    
    a_kin = np.loadtxt('et_%s_kin_en.xvg' % method)
    t_kin = a_kin[:,0]
    y_kin = a_kin[:,1]
    plt.plot(t_pot, y_pot, "bo", t_kin, y_kin, "ro", alpha=0.8)
    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('Energies (%s method)' % method) 
    plt.show()

Видно, что в большинстве случаев значения обеих энергий практически совпадают. Исключение составляет метод nh, где кинетическая энергия явно выше потенциальной. По форме распределения энергий метод be отличается от всех остальных: амплитуда колебаний значений уменьшается со временем. У остальных методов она остается примерно постоянной, а различаются только абсолютные значения амплитуды: самые маленькие у an, самые большие у nh.

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

Сначала создадим индекс файл с одной связью (в текстовом редакторе, файл назовём b.ndx):

[ b ] 1 2

Запустим утилиту по анализу связей g_bond:

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

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

In [42]:
for method in ["be","vr","nh","an","sd"]:
    a = np.loadtxt('bond_%s.xvg' % method)
    x = a[:,0]
    y = a[:,1]
    width = 0.0003
    plt.bar(x, y, width, color="green")
    plt.title('Length of C-C bond (%s method)' % method)
    plt.show()

Распределению Максвелла-Больцмана наиболее соответствуют распределения длин С-С связей в системах vr и sd. В системах nh и be большое количество молекул имеют длину связи, близкую к среднему значению, длина С-С связи меняется незначительно, что, вероятно, менее правдоподобно.

Выводы:

Метод Андерсена (an) в нашем случае производит только изменение длины С-С связи, имеет очень маленькую амплитуду значений энергий молекулы, для моделирования системы из одной молекулы не подходит. Метод Нуза-Хувера (nh) даёт сильные выбросы в значениях энергии (особенно кинетической) и не сильно меняет длину С-С связи. При использовании метода Берендсена (be) амплитуда изменения энергий постепенно снижается, что кажется не очень правдоподобным. По характерам изменения энергий и длин С-С связей наиболее правдоподобными кажутся метод стохастической молекулярной динамики (sd) и метод "Velocity rescale" (vr), но метод sd даёт очень хаотичные изменения в положении молекулы и помимо этого является самым медленным, поэтому, на мой взгляд, оптимальным в данном случае является метод "Velocity rescale" (vr).