#Получим файл с координатами для молекулы этана
! 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
Данные программы были запущены через putty на kodomo.
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)
Метод Берендсена показывает лучшую сходимость. Метод Андерсена демонстрирует минимальные колебания, но относитльно его сходимости трудно сделать вывод. Метод Нуза-Хувера демонстрирует наибольшую амплитуду колебаний энергии, трудно сделать вывод относительно его пригодности для рассчетов. Методы Velocity rescale и стохастической молекулярной динамики демонстрируют похожие результаты, но хуже, чем алгоритм Берендсена
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
Данные программы были запущены через putty на kodomo.
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)