Будем изучать как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования это одна молекула этана. Начнем с того, что подготовим файл координат и файл топологии.
import subprocess
topheader = '''#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_135 1 ETH C1 1 -0.189 12.01
2 opls_135 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
'''
bondheader='''
[ bonds ]
; ai aj funct b0 kb
'''
angheader='''
[ angles ]
; ai aj ak funct phi0 kphi
3 1 4 1
4 1 5 1
3 1 5 1
2 1 3 1
2 1 4 1
2 1 5 1
'''
diheader = '''
[ dihedrals ]
; ai aj ak al funct
3 1 2 6 3
3 1 2 7 3
3 1 2 8 3
'''
pairs='''
[ pairs ]
; 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
'''
top = open('./et.top', 'w')
top.write(topheader)
top.write(bondheader)
top.write(' 1 2 1'+'\n')
for i in range(3,6):
top.write(' 1 '+str(i)+' 1'+'\n')
for i in range(6,9):
top.write(' 2 '+str(i)+' 1'+'\n')
top.write('\n')
top.write(angheader)
for i in range(6,9):
top.write(' '+str(i)+' 2 1 1'+'\n')
for i in range(6,9):
for j in range(i+1,9):
top.write(' '+str(i)+' 2 '+str(j)+' 1'+'\n')
top.write('\n')
top.write(diheader)
for i in range(4,6):
for j in range(6,9):
top.write(' '+str(i)+' 1 2 '+str(j)+' 3'+'\n')
top.write('\n')
top.write(pairs)
top.close()
Запустим динамики молекулы этана с разными параметрами контроля температуры - V-rescale (vr), Berendsen (be), nose-hoover (nh), Andersen (an), без контроля температуры (sd).
def run_gmx(typea):
do = 'grompp -f '+typea+'.mdp -c et.gro -p et.top -o ./et_'+typea+'.tpr'
p = subprocess.Popen(do,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out=p.communicate()[0]
with open('grompp'+typea+'.out', 'w') as outfile:
outfile.write(out)
do2 = 'mdrun -deffnm et_'+typea+' -v -nt 1'
p1 = subprocess.Popen(do2,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out2=p1.communicate()[0]
with open('mdrun'+typea+'.out', 'w') as outfile2:
outfile2.write(out2)
list_of_ts = ['nh','sd','vr','be','an']
#for i in list_of_ts:
# run_gmx(i)
for i in list_of_ts:
do = 'echo 0 | trjconv -f et_'+i+'.trr -s et_'+i+'.tpr -o et_'+i+'.pdb'
p = subprocess.Popen(do,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
Посмотрев на полученные траектории можно заметить, что наиболее хаотично движется молекула этана из траектории без контроля температуры. Наиболее статична молекула в траектории an, не наблюдается вращения молекулы, присутствует только колебание длин связей. В траектории nh наблюдается вращение молекулы вкруг оси, проходящей через связь СС. В траекториях be и vr также есть вращение по другим осям в пространстве.
Построим графики изменения потенциальной и кинетической энергии молекулы в течение симуляции.
for i in list_of_ts:
do = 'echo 10 | g_energy -f et_'+i+'.edr -o et_'+i+'_pot.xvg -xvg none'
p = subprocess.Popen(do,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
do2 = 'echo 11 | g_energy -f et_'+i+'.edr -o et_'+i+'_kin.xvg -xvg none'
p = subprocess.Popen(do2,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
do2 = 'echo 12 | g_energy -f et_'+i+'.edr -o et_'+i+'_tot.xvg -xvg none'
p = subprocess.Popen(do2,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
import numpy as np
import matplotlib.pyplot as plt
def plot_ens(typea):
a= np.loadtxt('et_'+typea+'_pot.xvg')
t=a[:,0]
y=a[:,1]
plt.plot(t,y)
def plot_kin(typea):
a= np.loadtxt('et_'+typea+'_kin.xvg')
t=a[:,0]
y=a[:,1]
plt.plot(t,y)
plt.rcParams["figure.figsize"] = [15,3]
for i in list_of_ts:
plot_ens(i)
plt.legend(list_of_ts, loc='best')
plt.show()
plt.rcParams["figure.figsize"] = [15,3]
for i in list_of_ts:
plot_kin(i)
plt.legend(list_of_ts, loc='best')
plt.show()
plt.rcParams["figure.figsize"] = [15,3]
for i in list_of_ts:
plot_ens(i)
plt.legend(i+'pot')
axes = plt.gca()
axes.set_ylim([0,300])
plt.show()
plot_kin(i)
plt.legend(i+'kin')
axes = plt.gca()
axes.set_ylim([0,300])
plt.show()
Термостат Андерсена удерживает и кинетическую, и потенциальную энергию системы на минимуме на протяжении всей симуляции. Термостат Берендсена в течение симуляции гасит флуктуации энергии системы. Алгоритм V-rescale дает такое же распределение по энергии системы, как симуляция без контроля температуры. Термостат nose-hoover позволяет энегрии системы колебаться сильнее, чем вообще без контроля температуры.
Построим распределение длины связи СС в течение симуляции.
b = '''[ b ]
1 2 '''
with open('b.ndx', 'w') as outfilea:
outfilea.write(b)
for i in list_of_ts:
do = 'g_bond -f et_'+i+'.trr -s et_'+i+'.tpr -o bond_'+i+'.xvg -n b.ndx -xvg none'
p = subprocess.Popen(do,
shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
def plot_cc(typea):
a= np.loadtxt('bond_'+typea+'.xvg')
x=a[:,0]
n=a[:,1]
plt.plot(x,n)
for i in list_of_ts:
plot_cc(i)
axes = plt.gca()
axes.set_xlim([0.14,0.17])
plt.legend(i+'cc')
plt.show()
#plt.hist(cvs[i], bins, alpha=0.2, color = col, edgecolor=col)
Судя по форме распределения длины связи, термостат Берендсена справляется с задачей реалистичного поддержания температуры системы лучше всех.