Pr7: Вычисление параметров для молекулярной механики
In [ ]:
# Все нужные файлы были скачаны с сервера kodomo.
# Файл с топологией дополнили для дальнейшей работы. Тип для торсионных углов - 3.
In [ ]:
# В kodomo провели симуляции со всеми 5ю предложенными методами контроля температуры, затем получили pdb-файлы,
# которые тоже скачали
In [2]:
import numpy as np
import matplotlib.pyplot as plt
In [3]:
data = ['be', 'an', 'vr', 'nh', 'sd']
text = ['Berendsen', 'Anderson', 'Velocity-rescale', 'Nooze-Hoover', 'Stochastic MD']
In [4]:
# создаем файл с одной С-С связью
with open('b.ndx', 'w') as a:
a.write('''[ b ]
1 2 ''')
In [7]:
# Рисуем графики потенциальной и кинетической энергий
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='orange', alpha=0.9)
plt.plot(t, y_kin, c='yellow', alpha=0.6)
plt.title(txt)
plt.legend(['E (full)', 'E (potential)', 'E (kinetic)'])
plt.show()
In [8]:
# Рисуем графики распределения потенциальной и кинетической энергий
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='orange')
plt.title(txt)
plt.show()
# Берендсен и Нуз-Хувер больше всего напоминают Больцмановское распределение
In [9]:
# Рисуем распределение длины связи C-C
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='red')
plt.title(txt)
plt.show()
#Вышли вполне нормальные распределения.
#Длина С-С - в среднем 0,153 нм
In [11]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
In [12]:
import pymol
pymol.finish_launching()
from pymol import cmd
In [13]:
# посмотрим наконец на результаты моделей
for d in data:
cmd.load(f"et_{d}.pdb")
