Создадим файл топологии для молекулы этана. Для этого был взят шаблон с сайта кодомо, который был дополнен связями, углами и торсионными углами, также из файла oplsaa.ff/atomtypes.atp были взяты типы атомов для С и Н (у алканов), а из файла oplsaa.ff/forcefield.itp мы выяснили, что тип функционала для торсионных углов равен 3, а не 1.

%%bash

cat /usr/share/gromacs/top/oplsaa.ff/ffbonded.itp
cat /usr/share/gromacs/top/oplsaa.ff/ffbonded.itp

Получился такой файл et.top.

#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
    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

#метод Берендсена
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp
#метод "Velocity rescale"
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp
#метод Нуза-Хувера
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp
#метод Андерсена
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp
#метод стохастической молекулярной динамики
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp

Создаем входные файлы для молекулярной динамики.

%%bash
cd pr7_pymol
for i in *mdp
do
s=$(echo $i | cut -d '.' -f 1)
echo $s
grompp -f $i -c et.gro -p et.top -o et_$s.tpr
done

Запускаем саму программу mdrun.

%%bash
cd pr7_pymol
for i in *tpr
do
mdrun -deffnm $i -v -nt 1
done

Теперь конвертируем результаты программы в pdb файлы и посмотрим на них в Pymol (данную команду пришлось вбивать вручную, так как она требует допонлительного ввода номера группы от пользователя).

%%bash
cd pr7_pymol
for i in *tpr
do
trjconv -f $i.trr -s $i -o $i.pdb
done

Локально был запущен приведенный ниже скрипт (для каждого метода в отдельности).

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

import pymol
pymol.finish_launching()
from pymol import cmd,stored

cmd.do('''
delete all
load et_an.tpr.pdb
mset 1-251
set cache_frames=0
mclear
viewport 210,260
mpng et_an/mov
''')

import imageio
import os

folder = 'et_an' 
files = [f"{folder}\\{file}" for file in os.listdir(folder)]

images = [imageio.imread(file) for file in files]
imageio.mimwrite('et_an.gif', images, fps=24)

На выходе были получены следующие gif-анимации.

from IPython.display import Image
print('Метод Андерсена')
display(Image("pr7_pymol/et_an.gif", format='png'))
print('Метод Берендсена')
display(Image("pr7_pymol/et_be.gif", format='png'))
print('Метод Нуза-Хувера')
display(Image("pr7_pymol/et_nh.gif", format='png'))
print('Метод стохастической молекулярной динамики')
display(Image("pr7_pymol/et_sd.gif", format='png'))
print('Метод "Velocity rescale"')
display(Image("pr7_pymol/et_vr.gif", format='png'))

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

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

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

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

Метод "Velocity rescale"

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

Также для каждого метода были посчитаны потенциальная и кинетическая энергии (снова пришлось вбивать команду вручуную, так как она требует от пользователя указать номера параметров для выдачи, в нашем случае это потенциальная (10) и кинетическая (11) энергии).

%%bash
cd pr7_pymol
for i in *tpr
do
g_energy -f $i.edr -o $i_en.xvg -xvg none
done

import numpy as np
import matplotlib.pyplot as plt
print('Голубым обозначена потенциальная энергия, красным - кинетическая')
an= np.loadtxt('pr7_pymol/et_an_en.xvg')
plt.plot(an[:,0],an[:,1],'b-',an[:,0],an[:,2], 'r-')
plt.title('Метод Андерсена')
plt.xlabel('Время, пс')
plt.ylabel('Энергия, кДЖ/моль')
plt.show()
be= np.loadtxt('pr7_pymol/et_be_en.xvg')
plt.plot(be[:,0],be[:,1],'b-',be[:,0],be[:,2], 'r-')
plt.title('Метод Берендсена')
plt.xlabel('Время, пс')
plt.ylabel('Энергия, кДЖ/моль')
plt.show()
nh= np.loadtxt('pr7_pymol/et_nh_en.xvg')
plt.plot(nh[:,0],nh[:,1],'b-',nh[:,0],nh[:,2],'r-')
plt.title('Метод Нуза-Хувера')
plt.xlabel('Время, пс')
plt.ylabel('Энергия, кДЖ/моль')
plt.show()
sd= np.loadtxt('pr7_pymol/et_sd_en.xvg')
plt.plot(sd[:,0],sd[:,1],'b-',sd[:,0],sd[:,2],'r-')
plt.title('Метод стохастической молекулярной динамики')
plt.xlabel('Время, пс')
plt.ylabel('Энергия, кДЖ/моль')
plt.show()
vr= np.loadtxt('pr7_pymol/et_vr_en.xvg')
plt.plot(vr[:,0],vr[:,1],'b-',vr[:,0],vr[:,2],'r-')
plt.title('Метод "Velocity rescale"')
plt.xlabel('Время, пс')
plt.ylabel('Энергия, кДЖ/моль')
plt.show()

Голубым обозначена потенциальная энергия, красным - кинетическая

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

Теперь создадим индекс файл с одной связью (С-С) и посмотрим распределение ее длины в зависимости от метода.

%%bash
cd pr7_pymol
for i in *tpr
do
g_bond -f $i.trr -s $i -o bond_$i.xvg -n b.ndx -xvg none
done

an= np.loadtxt('pr7_pymol/bond_et_an.tpr.xvg')
plt.bar(an[:,0],an[:,1], width = 0.0005)
plt.title('Метод Андерсена')
plt.xlabel('Длина связи, нм')
plt.show()
be= np.loadtxt('pr7_pymol/bond_et_be.tpr.xvg')
plt.bar(be[:,0],be[:,1], width = 0.0005)
plt.title('Метод Берендсена')
plt.xlabel('Длина связи, нм')
plt.show()
nh= np.loadtxt('pr7_pymol/bond_et_nh.tpr.xvg')
plt.bar(nh[:,0],nh[:,1], width = 0.0005)
plt.title('Метод Нуза-Хувера')
plt.xlabel('Длина связи, нм')
plt.show()
sd= np.loadtxt('pr7_pymol/bond_et_sd.tpr.xvg')
plt.bar(sd[:,0],sd[:,1], width = 0.0005)
plt.title('Метод стохастической молекулярной динамики')
plt.xlabel('Длина связи, нм')
plt.show()
vr= np.loadtxt('pr7_pymol/bond_et_vr.tpr.xvg')
plt.bar(vr[:,0],vr[:,1], width = 0.0005)
plt.title('Метод "Velocity rescale"')
plt.xlabel('Длина связи, нм')
plt.show()

На графиках распределения длин связи С-С можно заметить, что у методов с меньшими колебаниями энергии также и меньший разброс значений этих длин (что в общем-то логично). При этом все методы дают в среднем длину связи 1.52-1.53 Ангстрема (что соответствует действительности).

Нарисуем распределения молекул по их потенциальным энергиям и сравим получившиеся распределения с распределением Больцмана.

an= np.loadtxt('pr7_pymol/et_an_en.xvg')
plt.hist(an[:,1], bins=100)
plt.title('Метод Андерсена')
plt.xlabel('Потенциальная энергия, кДж/моль')
plt.show()
be= np.loadtxt('pr7_pymol/et_be_en.xvg')
plt.hist(be[:,1], bins=100)
plt.title('Метод Берендсена')
plt.xlabel('Потенциальная энергия, кДж/моль')
plt.show()
nh= np.loadtxt('pr7_pymol/et_nh_en.xvg')
plt.hist(nh[:,1], bins=100)
plt.title('Метод Нуза-Хувера')
plt.xlabel('Потенциальная энергия, кДж/моль')
plt.show()
sd= np.loadtxt('pr7_pymol/et_sd_en.xvg')
plt.hist(sd[:,1], bins=100)
plt.title('Метод стохастической молекулярной динамики')
plt.xlabel('Потенциальная энергия, кДж/моль')
plt.show()
vr= np.loadtxt('pr7_pymol/et_vr_en.xvg')
plt.hist(vr[:,1], bins=100)
plt.title('Метод "Velocity rescale"')
plt.xlabel('Потенциальная энергия, кДж/моль')
plt.show()

Из выше представленных распределений видно, что метод Нуза-Хувера дает картину, больше всего напоминающую больцмановское распределение. На самом деле такой вид гистограмме придает тот факт, что при данном методе мы пробуем больший диапазон энергий (поэтому ось х идет вплоть до 250 кДЖ/моль). Тем не менее, большинство молекул тут действительно имеет низкие значения энергий (в отличии от других методов).

Сравним разные методы по количеству времени, затраченного на выполнение команды mdrun (real time, выдаваемое самим GROMACS).

import pandas as pd
data = [[3.913, 4.096, 4.085, 4.671, 4.147]]
pd.DataFrame(data, columns=["An", "Be", "Nh", "Sd", "Vr"])

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

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