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

Изучение реализации контроля температуры в молекулярной динамике на примере GROMACS; объект исследования - одна молекула этана.

Подготовим файл координат и файл топологии:

In [2]: 
%%bash
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/box_38.gro
make_ndx -f box_38.gro -o 1.ndx # выбираем остаток r1
editconf -f box_38.gro -o et1.gro -n 1.ndx # группа 3
#зададим ячейку и расположим молекулу по центру ячейку
editconf -f et1.gro -o et.gro -d 2 -c

Построим файл топологии 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
    1     2     6     1   
    6     2     8     1   
    6     2     7     1   
    7     2     8     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 ]
;  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

Подготовим входные файлы для mdrun:

In []: 
%%bash
m=(be vr nh an sd) # методы контроля температуры
for i in ${m[*]}; do grompp -f ${i}.mdp -c et.gro -p et.top -o et_${i}.tpr; done

Для каждого метода запустим mdrun:

In []: 
%%bash
for i in ${m[*]}; do mdrun -deffnm et_${i} -v -nt 1; done

Сконвертируем результаты в pdb и посмотрим в Pymol:

In []: 
%%bash
for i in ${m[*]}; do trjconv -f et_${i}.trr -s et_${i}.tpr -o et_${i}.pdb; done
In [3]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import time

cmd = ServerProxy(uri="http://localhost:9123/RPC2")

cmd.bg_color('white')
Dir = '/Users/aleksandrabezmenova/Documents/fbb/term8/Golovin/Pr7/'

ImPath = '/tmp/pymolpng.png'
def MakeImage():
    cmd.ray(500,400)
    cmd.png(ImPath)
    time.sleep(2)
    
def focus(sele):
    cmd.center(sele)
    cmd.zoom(sele)

Начальные конформации у всех пяти структур совпадают:

In [4]: 
MakeImage()
Image(ImPath)
Out[4]:

Ниже представолен структуры всех конформаций, полученных с помощью 5 методов:

In [5]: 
MakeImage()
Image(ImPath)
Out[5]:
In [6]: 
MakeImage()
Image(ImPath)
Out[6]:
In [7]: 
MakeImage()
Image(ImPath)
Out[7]:
In [8]: 
MakeImage()
Image(ImPath)
Out[8]:
In [9]: 
MakeImage()
Image(ImPath)
Out[9]:

an - молекула практически неподвижна

be - движется и вращается в ячейке, группы вращаются относительно С-С связи

nh - практически неподвижна в ячейке, группы вращаются относительно С-С связи

sd - движется хаотично

vr - движется и вращается в ячейке, группы вращаются относительно С-С связи

Неправдоподобным кажется результат an.

Сравним потенциальную энергию связи и кинетическую энергию для каждой из 5 систем:

In []: 
%%bash
for i in ${m[*]}; do g_energy -f et_${i}.edr -o et_${i}_en.xvg; done
In [44]: 
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [54]: 
def readFile(file):
    tt=[];uu=[];kk=[]
    for line in open(file,'r'):
        if line[0] in ["@","#"]: continue
        line=line.strip().split()
        tt.append(float(line[0]));uu.append(float(line[1]));kk.append(float(line[2]))
    return tt,uu,kk

fig, axes = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10,20))

methods=["an","be","nh","sd","vr"]

for i,m in enumerate(methods):
    tt,uu,kk=readFile("et_"+m+"_en.xvg")
    axes[i,0].plot(tt,uu,".")
    axes[i,1].plot(tt,uu,".")
    axes[i,0].set_xlabel("Time, ns")
    axes[i,0].set_title(m+"_U")
    axes[i,1].set_xlabel("Time, ns")
    axes[i,1].set_title(m+"_K")

plt.show()

an - энергия практически постоянна

be - энергия поначалу немного флуктуирует, затем флуктуации стихают, энергия приближается к 0

nh - в течении всего времени есть какое-то количество явных аутлайеров

sd, vr - постоянные не очень большие флуктуации вокруг значения ~20 кДж/моль

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

In [27]: 
%%bash
for i in ${m[*]}; do g_bond -f et_${i}.trr -s et_${i}.tpr -o bond_${i}.xvg -n b.ndx ; done
In [53]: 
def readFile(file):
    cc=[];pp=[]
    for line in open(file,'r'):
        if line[0] in ["@","#"]: continue
        line=line.strip().split()
        cc.append(float(line[0]));pp.append(float(line[1]))
    return cc,pp


fig, axes = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25,4))

for i,m in enumerate(methods):
    cc,pp=readFile("bond_"+m+".xvg")
    axes[i].bar(cc,pp,width=0.0001)
    axes[i].set_xlabel("C-C length, nm")
    axes[i].set_title(m)

plt.show()

Мне кажется, что наиболее реалистичные результаты дают методы Velocity rescale и стохастической молекулярной динамики (потенциальная и кинетическая энергия постоянно флуктуирует вогруг среднего; распределения длин связей похожи на больцмановское).