
#Получим файл с координатами для молекулы этана

! wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro

--2021-05-16 17:31:51--  http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 93.180.63.127
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:80... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro [following]
--2021-05-16 17:31:51--  https://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|93.180.63.127|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 399
Saving to: ‘etane.gro.1’

etane.gro.1         100%[===================>]     399  --.-KB/s    in 0s      

2021-05-16 17:31:51 (16.4 MB/s) - ‘etane.gro.1’ saved [399/399]


#Запишем файл топологии

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''')


programs = ['be', 'vr', 'nh', 'an', 'sd']


for i in programs:
    # скачаем конфигурационные файлы
    #!wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/{i}.mdp
    # строим входные файлы
    print('grompp -f '+i+'.mdp -c etane.gro -p et.top -o et_'+i+'.tpr')
    # запуск mdrun
    print('mdrun -deffnm et_'+i+' -v -nt 1')
    # конвертируем результаты в PDB
    print('trjconv -f et_'+i+'.trr -s et_'+i+'.tpr -o et_'+i+'.pdb')
    # сравним потенциальную и кинетическую энергии 
    print('g_energy -f et_'+i+'.edr -o et_'+i+'_en.xvg -xvg none')

grompp -f be.mdp -c etane.gro -p et.top -o et_be.tpr
mdrun -deffnm et_be -v -nt 1
trjconv -f et_be.trr -s et_be.tpr -o et_be.pdb
g_energy -f et_be.edr -o et_be_en.xvg -xvg none
grompp -f vr.mdp -c etane.gro -p et.top -o et_vr.tpr
mdrun -deffnm et_vr -v -nt 1
trjconv -f et_vr.trr -s et_vr.tpr -o et_vr.pdb
g_energy -f et_vr.edr -o et_vr_en.xvg -xvg none
grompp -f nh.mdp -c etane.gro -p et.top -o et_nh.tpr
mdrun -deffnm et_nh -v -nt 1
trjconv -f et_nh.trr -s et_nh.tpr -o et_nh.pdb
g_energy -f et_nh.edr -o et_nh_en.xvg -xvg none
grompp -f an.mdp -c etane.gro -p et.top -o et_an.tpr
mdrun -deffnm et_an -v -nt 1
trjconv -f et_an.trr -s et_an.tpr -o et_an.pdb
g_energy -f et_an.edr -o et_an_en.xvg -xvg none
grompp -f sd.mdp -c etane.gro -p et.top -o et_sd.tpr
mdrun -deffnm et_sd -v -nt 1
trjconv -f et_sd.trr -s et_sd.tpr -o et_sd.pdb
g_energy -f et_sd.edr -o et_sd_en.xvg -xvg none


import numpy as np
import matplotlib.pyplot as plt
import time


def plot_energies(program):
    en = np.loadtxt('et_%s_en.xvg' % program)
    x = en[:,0]
    potential = en[:,1]
    kinetic = en[:,2]
    plt.figure(figsize=(10, 6))
    plt.plot(x, potential, color='red', alpha=0.5)
    plt.plot(x, kinetic, color='blue', alpha=0.5)
    plt.ylim(0, 500)
    plt.legend(['Potential energy', 'Kinetic energy'])
    plt.title(program + " energies")


for program in programs:
    plot_energies(program)


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


for i in programs:
        # оцениваем длину C-C связи
    print('g_bond -f et_'+i+'.trr -s et_'+i+'.tpr -o bond_'+i+'.xvg -n b.ndx -xvg none')

g_bond -f et_be.trr -s et_be.tpr -o bond_be.xvg -n b.ndx -xvg none
g_bond -f et_vr.trr -s et_vr.tpr -o bond_vr.xvg -n b.ndx -xvg none
g_bond -f et_nh.trr -s et_nh.tpr -o bond_nh.xvg -n b.ndx -xvg none
g_bond -f et_an.trr -s et_an.tpr -o bond_an.xvg -n b.ndx -xvg none
g_bond -f et_sd.trr -s et_sd.tpr -o bond_sd.xvg -n b.ndx -xvg none


def bondlen_plot(program, w):
    tabs = np.loadtxt('bond_' + program + '.xvg')
    x = tabs[:,0]
    y = tabs[:,1]
    plt.figure(figsize=(10, 6))
    plt.bar(x, y, w, color='lightgreen')
    plt.title('C-C bond length distribution for the ' + program + ' algorithm')


def distr_cclen(program):
    cclen = np.loadtxt('bond_%s.xvg' % program)
    x = cclen[:,0]
    y = cclen[:,1]
    plt.figure(figsize=(10, 6))
    plt.bar(x, y, 0.0005, color='lightgreen')
    plt.title(program + ' C-C bond lengths')


for program in programs:
    distr_cclen(program)

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

Сравнение алгоритмов контроля температур¶