Практикум 7. Вычисление параметров для молекулярной механики

Скопировали в рабочую папку все исходные файлы с сервера kodomo.

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

Сам я думал изначально на другие типы opls, но с ними ничего не работало. Исправил на эти (139, 140).

Тип для торсионных углов 3.

In [3]:
%cat 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_139      1    ETH      C1      1    -0.189      12.011
    2   opls_139      1    ETH      C2      2    -0.155      12.011
    3   opls_140      1    ETH      H1      3     0.060       1.008
    4   opls_140      1    ETH      H2      4     0.060       1.008
    5   opls_140      1    ETH      H3      5     0.060       1.008
    6   opls_140      1    ETH      H4      6     0.060       1.008
    7   opls_140      1    ETH      H5      7     0.060       1.008
    8   opls_140      1    ETH      H6      8     0.060       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

Посчитал работу симуляции при разном контроле температуры (вручную)

Сначала:

  1. Строил входные файлы для mdrun с помощью grompp.
  2. Запускал на каждом из файлов mdrun.
  3. Каждый из результатов сконвертировал в pdb (trjconv).
  4. И сравнил потенциальную и кинетическую энергию (g_energy).
  5. Также создал файл b.ndx ([ b ]\n1 2) для выдачи распределений длины С-С связи за моделирование (g_bond).

Сделав все это, перенес данные на свой компьютер.

In [11]:
import numpy as np
import matplotlib.pyplot as plt
In [13]:
data = ['be', 'an', 'vr', 'nh', 'sd']
text = ['Berendsen', 'Anderson', 'Velocity-rescale', 'Nooze-Hoover', 'Stochastic MD']

Нарисуем энерг. сост. (для потенциальной и кинетической)

In [25]:
plt.subplots(2, 3, figsize=(16, 8))
for i, (d, txt) in enumerate(zip(data, text), 1):
    energy = np.loadtxt(f'et_{d}_en.xvg')
    t = energy[:,0]
    y_pot = energy[:,1]
    y_kin = energy[:,2]
    
    plt.subplot(2, 3, i)
    
    plt.plot(t, y_kin + y_pot, c='gray', alpha=0.2)
    plt.plot(t, y_pot, c='green', alpha=0.9)
    plt.plot(t, y_kin, c='chartreuse', alpha=0.6)
    plt.title(txt)
    plt.legend(['E (full)', 'E (potential)', 'E (kinetic)'])
plt.show()

Берендсен сходится лучше всех.

Распределение энергии:

In [37]:
plt.subplots(2, 3, figsize=(16, 8))
for i, (d, txt) in enumerate(zip(data, text), 1):
    energy = np.loadtxt(f'et_{d}_en.xvg')
    e_pot = energy[:,1]
    
    plt.subplot(2, 3, i)
    
    plt.hist(e_pot, bins=100, color='violet')
    plt.title(txt)
plt.show()

Мне кажется, на распределение Максвелла больше всего похож стохастический метод и, возможно, V-rescale.

И распределение длины связи (1-2) по времени

In [32]:
plt.subplots(2, 3, figsize=(16, 8))
for i, (d, txt) in enumerate(zip(data, text), 1):
    bond = np.loadtxt(f'et_{d}_bond.xvg')
    x = bond[:,0]
    y = bond[:,1]
    
    plt.subplot(2, 3, i)
    
    plt.bar(x, y, width=0.0002, color='salmon')
    plt.title(txt)
plt.show()

Моих мозгов вначале хватило, чтобы сказать себе: "Хы, на ГЗ МГУ похоже".

Ну и посмотрим на все это в pymol

In [38]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
In [39]:
import pymol
pymol.finish_launching()
from pymol import cmd
In [40]:
for d in data:
    cmd.load(f"et_{d}.pdb")

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

В стохастическом же методе молекула смещает свой центр масс без внешней энергии, это тоже, на мой взгляд, нехорошо

В методе Андерсона молекула фактически не меняла конформацию, а вот в методе V-R были некоторые переклики из одной анти-конформации в другую, что больше соответствует реальному поведению.

Такой визуальный анализ, конечно, не несет глубокого смысла, особенно с учетом того, что я мог задать не те параметры на входе =(

Вот, кажется, и все. Я бы предпочел термостат Velocity-rescale остальным, но недаром их столько сделали - применение зависит от задачи.