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

Необходимо получить файл с топологией этана. Для этого нужно взять заготовку файла et.top из данного задания и заменить типы атомов и значения углов связей в соответствии с данными из файла atomtypes.atp на kodomo.

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
import os
import subprocess
import matplotlib.pyplot as plt

# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )

0

# и придумаем циклы 
## связи
print('Bonds')
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
    print('    ',b.GetBeginAtomIdx()+1 ,'    ',b.GetEndAtomIdx()+1,'    1')
## углы 
print('Angles')
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')
## dihedrals    
print('Dihedrals')
for b1 in m3d.GetBonds():
    for b2 in m3d.GetBonds():
        for b3 in m3d.GetBonds():
            if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx() and b2.GetEndAtomIdx()==b3.GetBeginAtomIdx() and b3.GetIdx() != b2.GetIdx():
                print(' ',b1.GetEndAtomIdx()+1,'   ',b1.GetBeginAtomIdx()+1,'   ',b2.GetEndAtomIdx()+1,'   ',b3.GetEndAtomIdx()+1,'   ','3')

Bonds
     1      2     1
     1      3     1
     1      4     1
     1      5     1
     2      6     1
     2      7     1
     2      8     1
Angles
     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
Dihedrals
  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

В качестве типов атомов были взяты opls_139 (alkane C) и opls_140 (alkane H)

#запишем результаты в файл топологии
with open('et.top', 'w') as f:
    f.write('''#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_139      1    ETH      C1      1    -0.189      12.01
    2   opls_139      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''')

Даны 5 конфигурационных файлов с разными алгоритмами контроля температуры:
be.mdp - метод Берендсена для контроля температуры.
vr.mdp - метод "Velocity rescale" для контроля температуры.
nh.mdp - метод Нуза-Хувера для контроля температуры.
an.mdp - метод Андерсена для контроля температуры.
sd.mdp - метод стохастической молекулярной динамики.
Все файлы поместила с рабочую директорию.

for program in ['be','vr','nh','an','sd']:
    #построение входных файлов
    cmd = 'grompp -f %s.mdp -c ethane.gro -p et.top -o et_%s.tpr >& log.grompp' %(program,program)
    subprocess.call(cmd,shell=True)
    #запуск молекулярно-динамического движка mdrun
    cmd = 'mdrun -deffnm et_%s.tpr -v -nt 1' %program
    subprocess.call(cmd,shell=True)
    #конвертация в pdb
    cmd = 'echo 0 | trjconv -f et_%s.tpr.trr -s et_%s.tpr -o et_%s.pdb >& log.trjconv' %(program,program,program)
    subprocess.call(cmd,shell=True)

Посмотрим на итоговые pdb в PyMol

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

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

cmd.do('''
delete all
load et_an.pdb
mset 1-250
set cache_frames=0
mclear
viewport 300,300
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.mimsave('et_an.gif', images, fps=24)

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

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

'Mетод "Velocity rescale"'

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

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

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

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

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

import subprocess
for program in ['be','vr','nh','an','sd']:
    cmd1 = 'echo 10 | g_energy -f et_%s.tpr.edr -o et_%s_pot_en.xvg -xvg none >& log.g_energy' %(program,program)
    cmd2 = 'echo 11 | g_energy -f et_%s.tpr.edr -o et_%s_kin_en.xvg -xvg none >& log.g_energy' %(program,program)
    subprocess.call(cmd1,shell=True)
    subprocess.call(cmd2,shell=True)

data = ['be', 'an', 'vr', 'nh', 'sd']
text = ['Berendsen', 'Anderson', 'Velocity-rescale', 'Nooze-Hoover', 'Stochastic molecular dinamics']

tex = 0
fig,axs = plt.subplots(2, 3, figsize=(16, 8))
for d in data:
    a= np.loadtxt('et_%s_pot_en.xvg' % d)
    b= np.loadtxt('et_%s_kin_en.xvg' % d)
    t1=a[:,0]
    y1=a[:,1]
    t2=b[:,0]
    y2=b[:,1]
    print([tex // 3, tex%3])
    axs[tex // 3, tex%3].plot(t1, y1, "blue", label = 'Potential energy', alpha=0.4)
    axs[tex // 3, tex%3].plot(t2, y2, "pink", alpha=0.4, label = 'Kinetic energy')
    axs[tex // 3, tex%3].set_title('Energies of %s system' % text[tex]) 
    axs[tex // 3, tex%3].legend()
    axs[tex // 3, tex%3].set_xlabel('t, sec')
    tex +=1
plt.tight_layout()

[0, 0]
[0, 1]
[0, 2]
[1, 0]
[1, 1]

В случае методов Берендсена и Нуза-Хувера наблюдаются колебания с уменьшающейся амлитудой (у nh размах намного больше и со временем молекула очень теряет в энергии). Метод Андерсена дает постоянные по силе колебания с очень маленькой амплитудой энергии - на gif молекула почти не двигалась. Картины изменения энергий в стохастическом и velocity rescale случаях очень похожи.

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

with open('b.ndx', 'w') as a:
    a.write('''[ b ]
1 2 ''')

import subprocess
for program in ['be','vr','nh','an','sd']:
    cmd = 'g_bond -f et_%s.tpr.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none' % (program,program,program)
    subprocess.call(cmd,shell=True)

tex = 0
fig,axs = plt.subplots(2, 3, figsize=(16, 8))
for d in data:
    bond= np.loadtxt('bond_%s.xvg' % d)
    x=bond[:,0]
    y=bond[:,1]
    axs[tex // 3, tex%3].bar(x, y, width=0.0003)
    axs[tex // 3, tex%3].set_title(text[tex]) 
    #axs[tex // 3, tex%3].legend()
    axs[tex // 3, tex%3].set_xlabel('Bond length, nm')
    tex +=1
fig.suptitle('С-С bond length distribution', fontsize=16, x = 0.5, y = 1.1)
plt.tight_layout()

Размах значений по длине связей мал у метово с небольшим разбросом энергии.
Длина С-С связи равна 0.1522 нм, около этого значения и колеблются максимумы всех методов.

Посмотрим распределение по потенциальным энергиям

tex = 0
fig,axs = plt.subplots(2, 3, figsize=(16, 8))
for d in data:
    pot= np.loadtxt('et_%s_pot_en.xvg' % d)
    x=pot[:,0]
    y=pot[:,1]
    axs[tex // 3, tex%3].hist(y, 100)
    axs[tex // 3, tex%3].set_title(text[tex]) 
    tex +=1
fig.suptitle('Potential energy distribution', fontsize=16, x = 0.5, y = 1.1)
plt.tight_layout()

На распределение Больцмана больше всего похож график метода Нуза-Хувера.

Сравним время работы алгоритмов

tex = 0
print ('               NODE (s)   Real (s)      (%)')
for item in data:
    with open('et_%s.tpr.log' % item, 'r') as log:
        handle = log.readlines()
        for line in range(len(handle)):
            if 'NODE (s)   Real (s)      (%)' in handle[line]:
                print (text[tex])
                print(handle[line+1])
    tex +=1

               NODE (s)   Real (s)      (%)
Berendsen
       Time:      3.870      4.033     96.0

Anderson
       Time:      3.750      4.044     92.7

Velocity-rescale
       Time:      3.900      4.084     95.5

Nooze-Hoover
       Time:      3.940      4.083     96.5

Stochastic molecular dinamics
       Time:      4.410      4.691     94.0

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