Запустим PyMol для работы в режиме сервера, а также определим несколько функций.

from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys, time

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

# __main__.pymol_argv = [ 'pymol', '-cp' ]  # to disable GUI
import pymol
from pymol import cmd
pymol.finish_launching()
 
from IPython.display import Image

defaultImage = '/tmp/pymolimg.png'
def prepareImage(width=300, height=300, sleep=2, filename=defaultImage):
    ## To save the rendered image
    cmd.ray(width, height)
    cmd.png('/tmp/pymolimg.png')
    time.sleep(sleep)

def showPDBoneState(name):
    cmd.do('''
    reinit
    bg white
    set antialias, 2
    set ray_trace_mode, 3
    load ''' + name + '.pdb' + '''
    split_states ''' + name + '''
    hide everything, all
    select only, ''' + name + '_0001' + '''
    show spheres, only
    set sphere_scale, 0.2
    show sticks, only
    set stick_radius, 0.1
    zoom only
    ''')

    from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys, time

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

# __main__.pymol_argv = [ 'pymol', '-cp' ]  # to disable GUI
import pymol
from pymol import cmd
pymol.finish_launching()
 
from IPython.display import Image

defaultImage = '/tmp/pymolimg.png'
def prepareImage(width=300, height=300, sleep=2, filename=defaultImage):
    ## To save the rendered image
    cmd.ray(width, height)
    cmd.png('/tmp/pymolimg.png')
    time.sleep(sleep)

# Define some shortcuts

def showPDBoneState(name):
    cmd.do('''
    reinit
    bg white
    set antialias, 2
    set ray_trace_mode, 3
    load ''' + name + '.pdb' + '''
    split_states ''' + name + '''
    hide everything, all
    select only, ''' + name + '_0001' + '''
    show spheres, only
    set sphere_scale, 0.2
    show sticks, only
    set stick_radius, 0.1
    zoom only
    ''')

def createPNGsFromStates(name):
    cmd.do('''
    reinit
    bg white
    set antialias, 2
    set ray_trace_mode, 3
    load ''' + name + '.pdb' + '''
    hide everything, all
    show spheres, all
    set sphere_scale, 0.2
    show sticks, all
    set stick_radius, 0.1
    color iridium, e. C
    color teal, e. H''')
    nstates = cmd.count_states("(all)")
    cmd.mset('1-'+str(nstates))
    cmd.mpng(name+'_state')

from IPython.display import HTML
from base64 import b64encode
def showMovie(filename):
    video = open(filename, 'rb').read()
    video_encoded = b64encode(video)
    video_tag = '<video controls alt="PyMol Movie" src="data:video/mp4;base64,{0}" type="video/mp4">'.format(video_encoded)
    return HTML(data=video_tag)

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

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

Сначала подготовим файл координат и файл топологии. На основе .gro файла с 38 молекулами этана создадим индекс-файл, в котором будет группа из одной молекулы этана.

%%bash
make_ndx -f box_38.gro -o 1.ndx
# Then type `r1` to choose the first residue
# Then type `q` to save and quit

Теперь создадим .gro файл с одной молекулой и зададим ячейку.

При запуске ediconf выберем номер соответствующей группы из одной молекулы (3). Затем зададим ячейку и расположим молекулу по центру ячейки.

%%bash
editconf -f box_38.gro -o et1.gro -n 1.ndx
editconf -f et1.gro -o et.gro -d 2 -c

Построим файл топологии et.top для этана. При этом необходимо, в частности, подобрать типы атомов.

ethane = """#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   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      1  
    3    1     2     7      1 
    3    1     2     8      1  
    4    1     2     6      1  
    4    1     2     7      1  
    4    1     2     8      1  
    5    1     2     6      1 
    5    1     2     7      1  
    5    1     2     8      1  
[ 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"""

with open('et.top', 'w') as w:
    w.write(ethane)

Далее будем работать с пятью файлами с разными параметрами контроля температуры:

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

Сначала построим входные файлы для молекулярно-динамического движка mdrun с помощью grompp, а затем запустим сам движок mdrun.

%%bash
for i in $(echo be vr nh an sd); do
    grompp -f $i.mdp -c et.gro -p et.top -o et_$i.tpr
    mdrun -deffnm et_$i -v -nt 1;
done

Анализ результатов

Проанализируем результаты. Сначала с помощью PyMol визуализируем каждую из систем.

%%bash
for i in $(echo be vr nh an sd); do
    trjconv -f et_$i.trr -s et_$i.tpr -o et_$i.pdb;
done

Геометрические характеристики начальных состояний систем идентичны.

showPDBoneState('et_be')
# showPDBoneState('et_be')
# showPDBoneState('et_vr')
# showPDBoneState('et_nh')
# showPDBoneState('et_an')
# showPDBoneState('et_sd')
prepareImage(); Image(defaultImage)

Различия выявляются в последовательности состояний (states) в PyMol, соответствующих состояниям системы в разные моменты времени.

Наиболее отлично от других применение метода Нуза-Хувера: в системе изменяется только длина одной C — H связи. При использовании других методов изменяются как длины всех связей, так и углы молекулы.

createPNGsFromStates('et_vr')

%%bash
ffmpeg -i et_vr_state0%03d.png -c:v libx264 -pix_fmt yuv420p \
    -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2,setpts=3.0*PTS" et_vr_states.mp4 2> /dev/null

Например, в случае использования метода Velocity rescale это выглядит следующим образом.

showMovie('et_vr_states.mp4')

Сравним потенциальную энергию связи и кинетическую энергию для каждой из 5 систем. Для этого после запуска команды будем выбирать параметры под номерами 10 и 11.

%%bash
for i in $(echo be vr nh an sd); do
    g_energy -f et_$i.edr -o et_$i_en.xvg;
done

Построим графики изменения энергий.

%matplotlib inline
import matplotlib.pyplot as plt
import math

def parseXVG(filename):
    with open(filename, 'r') as f:
        lines = filter(lambda x: not (x.startswith('#') or x.startswith('@')), f.readlines())
        [t, V, T] = zip(*[map(float, line.split()) for line in lines])
        return [t, V, T]  # [time (ps), potential energy (kJ/mol), kinetic energy (kJ/mol)]

def liftF(f, xs, ys):
    return f(map(f, [xs, ys]))

def log10vec(vector):
    return map(lambda x: math.log(x, 10), vector)

def showEnergy(filename, descr="", lpos="lower right"):
    [x, y, z] = parseXVG(filename)
    
    f = plt.figure(1)
    f.text(0.5, 1, descr, 
           horizontalalignment='center', verticalalignment='top')
    
    plt.subplot(141)
    plt.title("Potential energy")
    plt.scatter(x, y, c='b', alpha=0.7)
    plt.xlabel("Time, ps")
    plt.ylabel("Energy, kJ/mol")
    
    plt.subplot(142)
    plt.title("Kinetic energy")
    plt.scatter(x, z, c='r', alpha=0.7)
    plt.legend(loc=lpos)
    plt.xlabel("Time, ps")
    
    plt.subplot(143)
    plt.title("log10 potential energy")
    plt.scatter(x, log10vec(y), c='c', alpha=0.7)
    plt.legend(loc=lpos)
    plt.xlabel("Time, ps")
    plt.ylabel("log10 energy, kJ/mol")
    
    plt.subplot(144)
    plt.title("log10 kinetic energy")
    plt.scatter(x, log10vec(z), c='m', alpha=0.7)
    plt.legend(loc=lpos)
    plt.xlabel("Time, ps")
    
def showOneEnergy(listOfCoords, descr="", log=False, sparse=False):
    labels = ('be', 'vr', 'nh', 'an', 'sd')
    colours = ('r', 'g', 'b', 'y', 'c')
    f = plt.figure(1)
    f.text(0.5, 1, descr, 
       horizontalalignment='center', verticalalignment='top')

    for i, pairedList in enumerate(listOfCoords):
        (x, y) = pairedList
        plt.subplot(231+(i if i<4 else 5) if sparse else 151+i)
        plt.title(labels[i])
        if log:
            plt.scatter(x, log10vec(y), c=colours[i], alpha=0.7)
            if i==0: plt.ylabel("log10 Energy")
        else:
            plt.scatter(x, y, c=colours[i], alpha=0.7)
            if i==0 or i>2: plt.ylabel("Energy, kJ/mol")
    plt.xlabel("Time, ps")
    # plt.ylim(-.5*1e19, .5*1e19)
    plt.show()


filenames = ('et_be_en.xvg', 'et_vr_en.xvg', 'et_nh_en.xvg', 'et_an_en.xvg', 'et_sd_en.xvg')
coords = map(parseXVG, filenames)
getKinetic   = lambda k: (k[0], k[1])
getPotential = lambda k: (k[0], k[2])

Визуализируем изменение потенциальной энергии.

from matplotlib import pylab
pylab.rcParams['figure.figsize'] = (13.0, 5.0)
showOneEnergy(map(getPotential, coords), descr="Potential energy, kJ/mol")

Для наглядности лучше изобразим логарифм значения энергии.

showOneEnergy(map(getPotential, coords), descr="Potential energy, kJ/mol", log=True)

pylab.rcParams['figure.figsize'] = (15.0, 10.0)
showOneEnergy(map(getKinetic, coords), descr="Kinetic energy, kJ/mol", sparse=True)

Теперь построим графики изменений значений энергии для каждого метода отдельно.

pylab.rcParams['figure.figsize'] = (15.0, 4.0)
showEnergy('et_be_en.xvg', 'Berendsen`s method')

showEnergy('et_vr_en.xvg', 'Velocity rescale')

showEnergy('et_nh_en.xvg', 'Nose-Hoover method')

showEnergy('et_an_en.xvg', 'Andersen method')

showEnergy('et_sd_en.xvg', 'Stochastic method')

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

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

%%bash
echo -e "[ b ]\n1 2" > b.ndx
cat b.ndx

[ b ]
1 2

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

%%bash
for i in $(echo be vr nh an sd); do
    g_bond -f et_$i.trr -s et_$i.tpr -o bond_$i.xvg -d distance_$i.xvg -n b.ndx &> /dev/null;
done

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

labels = ('be', 'vr', 'nh', 'an', 'sd')
filenames = map(lambda x: 'distance_%s.xvg' % x, labels)
colours = ('r', 'g', 'b', 'y', 'c')


def parseXVG2d(filename):
    with open(filename, 'r') as f:
        lines = filter(lambda x: not (x.startswith('#') or x.startswith('@')), f.readlines())
        [t, d] = zip(*[map(float, line.split()) for line in lines])
        return [t, d]  # [time (ps), potential energy (kJ/mol), kinetic energy (kJ/mol)]

def distPlot(coords, type="scatter"):
    f = plt.figure(1)
    f.text(0.5, 1, "Distance", 
       horizontalalignment='center', verticalalignment='top')

    for i, pairedList in enumerate(coords):
        (x, y) = pairedList
        plt.subplot(231+(i if i<4 else 5))
        plt.title(labels[i])
        if type=="hist":
            plt.hist(y, color=colours[i], alpha=0.7)
        else:
            plt.scatter(x, y, color=colours[i], alpha=0.7)
        plt.xlabel("Time, ps")
        if i==0 or i>2: plt.ylabel("Distance, nm")
    plt.xlabel("Time, ps")

distCoords = map(parseXVG2d, filenames)
distPlot(distCoords, "scatter")
pylab.rcParams['figure.figsize'] = (15.0, 10.0)

distPlot(distCoords, "hist")

Если исходить из всех наблюдений, то можно заключить, что наиболее реалистично поддерживать температуру в системе позволяют методы Velocity rescale и стохастической молекулярной динамики.

Учебный сайт

Бредихина Данилы

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

Анализ результатов

Ссылки на файлы