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

Данный практикум посвящен реализации контроля температуры в молекулярной динамике на примере GROMACS.

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

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

Объект исследования - одна молекула этана. Нам был предоставлен файл с координатами этана - etane.gro. Был построен файл с топологией et.top, используя описанный ниже скрипт:

In [1]:
# загрузим модули
import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem import Descriptors
from rdkit import RDConfig
from IPython.display import Image,display
import numpy as np
In [2]:
#создаем этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
Out[2]:
0
In [3]:
#прописываем связи между атомами
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
    print b.GetBeginAtomIdx()+1 , b.GetEndAtomIdx()+1, 1

1 2 1
1 3 1
1 4 1
1 5 1
2 6 1
2 7 1
2 8 1
In [4]:
#углы    
for b1 in m3d.GetBonds():
    for b2 in m3d.GetBonds():
        if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
            print b1.GetEndAtomIdx()+1 , b1.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1, 1

2 1 3 1
2 1 4 1
2 1 5 1
3 1 2 1
3 1 4 1
3 1 5 1
4 1 2 1
4 1 3 1
4 1 5 1
5 1 2 1
5 1 3 1
5 1 4 1
6 2 7 1
6 2 8 1
7 2 6 1
7 2 8 1
8 2 6 1
8 2 7 1
In [5]:
#торсионные углы
for b1 in m3d.GetBonds():
    for b2 in m3d.GetBonds():
        for b3 in m3d.GetBonds():
            if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b2.GetEndAtomIdx() == b3.GetBeginAtomIdx()and b1.GetIdx() != b2.GetIdx():
                print b1.GetEndAtomIdx()+1 , b1.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1, b3.GetEndAtomIdx()+1, 1

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

На основе файла ffbonded.itp, были выбраны правильные типы атомов: opls_139 для углерода и opls_140 для водорода.

In [6]:
import subprocess
import time
for method in ["be","vr","nh","an","sd"]:
    #С помощью grompp строим файлы .tpr для запуска молдинамики
    command1 = 'grompp -f %s.mdp -c etane.gro -p et.top -o et_%s.tpr' %(method,method)
    subprocess.Popen(command1,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    #запускаем молекулярную динамику
    start = time.time()
    command2 = 'mdrun -deffnm et_%s -v -nt 1' %method
    subprocess.Popen(command2,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    print("Runtime for %s-method: %s seconds" % (method,round(time.time() - start, 3)))

Runtime for be-method: 0.013 seconds
Runtime for vr-method: 0.013 seconds
Runtime for nh-method: 0.012 seconds
Runtime for an-method: 0.035 seconds
Runtime for sd-method: 0.026 seconds

Самыми быстрыми оказались методы be, vr, nh ~0,012-0,013 секунд, самыми медленными an и sd ~0,026-0,035 секунд.

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

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

In [7]:
#Переводим в файл с траекторий в .pdb формат и визуализируем:
for method in ["be","vr","nh","an","sd"]:
    command3 = 'echo 2 | trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb' %(method,method,method)
    subprocess.Popen(command3,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
In [31]:
Image(url='an.gif') #незначительные колебания относительно С-С связи, молекула не вращается
Out[31]:
In [26]:
Image(url='be.gif') #появляются относительно медленные вращетельные движения
Out[26]:
In [27]:
Image(url='nh.gif') #ускорение вращательного движения, подвижные CCH углы
Out[27]:
In [28]:
Image(url='sd.gif') #самая быстроизменяющаяся система, поступательные и вращательные движения
Out[28]:
In [29]:
Image(url='vr.gif') #поведение молекулы схоже с молекулой в системе be
Out[29]:

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

In [14]:
for method in ["be","vr","nh","an","sd"]:
    #создаем файл с потенциальной энергией системы
    command4 = 'echo 10 | g_energy -f et_%s.edr -o et_%s_en_potential.xvg -xvg none' %(method,method)
    subprocess.Popen(command4,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    #с кинетической
    command5 = 'echo 11 | g_energy -f et_%s.edr -o et_%s_en_kinetic.xvg -xvg none' %(method,method)
    subprocess.Popen(command5,shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
In [15]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
#для каждого метода сравниваем потенциальную и кинетическую энергии
for method in ["be","vr","nh","an","sd"]:
    a_potential= np.loadtxt('et_%s_en_potential.xvg' %method)
    a_kinetic= np.loadtxt('et_%s_en_kinetic.xvg' %method)
    t=a_kinetic[:,0]
    e_potential=a_potential[:,1]
    e_kinetic=a_kinetic[:,1]
    plt.plot(t, e_potential, color='grey')
    plt.plot(t, e_kinetic, color='orange')
    plt.title('Method %s' %method)
    kin_patch = mpatches.Patch(color='orange', label='Kinetic energy')
    pot_patch = mpatches.Patch(color='grey', label='Potential energy')
    plt.legend(handles=[kin_patch,pot_patch])
    plt.show()

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

Рассмотрим распределение длинны связи С-С за время моделирования. Сначала создадим индекс файл с одной связью. В текстовом редакторе создайте файл b.ndx со следующим содержимым:

[ b ]

1 2

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

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

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

In [24]:
for method in ["be","vr","nh","an","sd"]:
    bonds = np.loadtxt('bond_%s.xvg' % method)
    x = bonds[:,0]
    y = bonds[:,1]
    width = 0.0008
    plt.bar(x, y, width, color="grey")
    plt.title('Length of C-C bond (%s method)' % method)
    plt.xlim(0.14, 0.165)
    plt.show()

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

Вывод:

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