Переделанный файл:

for line in open('et.top', encoding='utf-8').readlines():
    print(line.strip())

#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

Сравнение молекулярной динамики¶

Я запустила все bash-скрипты. С помощью trjconv конвертировала результаты в pdb и посмотрела в PyMol. Динамика у всех моделей была разная и примерно может быть описана так:

Berendsen - 'Sympathy For The Devil', Rolling Stones
Velocity-rescale - 'Superfreak', Rick James
Nooze-Hoover - 'Heroin', Niko
Anderson - 'Нам с тобой', Виктор Цой
Stochastic MD - 'Rockit', Herbie Hancock

Сравнение потенциальной и кинетической энергии¶

import numpy as np
import matplotlib.pyplot as plt

algs = {'be': 'Berendsen', 'vr': 'Velocity-rescale', 'nh':'Nooze-Hoover', 'an':'Anderson', 'sd':'Stochastic MD'}

fig, axs = plt.subplots(5, figsize=(10, 20))
for i, d in enumerate(algs.keys()):
    energy = np.loadtxt(f'et_{d}.en.xvg')
    t = energy[:,0]
    e_pot = energy[:,1]
    e_kin = energy[:,2]
    axs[i].plot(t, e_pot, c='teal', alpha=0.7, label='Potential energy')
    axs[i].plot(t, e_kin, c='indianred', alpha=0.7, label='Kinetic energy')
    axs[i].legend()
    axs[i].set_title(algs[d])
plt.show()

В методах Velocity-rescale, Anderson, Stochastic MD обе энергии флуктуируют, а в Berendsen и Nooze-Hoover стабилизируются.

Сравнение длин связей за время моделирования¶

fig, axs = plt.subplots(5, figsize=(10, 18))
for i, d in enumerate(algs.keys()):
    x = np.loadtxt(f'bond_{d}.xvg')
    axs[i].bar(x[:,0], x[:,1], width=0.0005, color='teal')
    axs[i].set_title(algs[d])
    axs[i].set_xlim(0.145, 0.1625)
plt.show()

Распределение для Anderson какое-то вообще не реалистичное, у Stochastic MD и Velocity-rescale, наверное, дисперсия великовата, Berendsen и Nooze-Hoover, пожалуй, самые реалистичные.

Сравнение распределений по энергии¶

fig, axs = plt.subplots(5, figsize=(10, 18))
for i, d in enumerate(algs.keys()):
    energy = np.loadtxt(f'et_{d}.en.xvg')
    t = energy[:,0]
    e_pot = energy[:,1]
    e_kin = energy[:,2]
    axs[i].hist(e_pot + e_kin, bins=100, color='indianred')
    axs[i].set_title(algs[d])
plt.show()

На распределение Больцмана больше всего похоже распределение Nooze-Hoover. Учитывая все показатели, я бы назвала распределение Nooze-Hoover самым реалистичным.