Изучение работы методов контроля температуры в GROMACS

Начало задания был выполнено с помощью указанных ниже команд в консоли, так как так удобнее "общаться" с программами

In []: 
make_ndx -f box_38.gro -o 1.ndx
editconf -f box_38.gro -o et1.gro -n 1.ndx
editconf -f et1.gro -o et.gro -d 2 -c

Построили файл топологии et.top для этана, изменив типы атомов согласно /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp

Создадим скрипт, чтобы получить 5 файлов с разными параметрами контроля температуры:

◦be.mdp - метод Берендсена для контроля температуры;

◦vr.mdp - метод "Velocity rescale" для контроля температуры;

◦nh.mdp - метод Нуза-Хувера для контроля температуры;

◦an.mdp - метод Андерсена для контроля температуры;

◦sd.mdp - метод стохастической молекулярной динамики.

In []: 
%%bash
for i in {'be','vr','nh','an','sd'}; do 
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/${i}.mdp
grompp -f ${i}.mdp -c et.gro -p et.top -o et_${i}.tpr
mdrun -deffnm et_${i} -v -nt 1
echo "2" | trjconv -f et_${i}.trr -s et_${i}.tpr -o et_${i}.pdb
done

Результат в pdb получился. Все молекулы этана действительно молекулы этана. С разными параметрами контроля температуры по-разному выглядят возможные положения молекул (разрешены те или иные степени свободы)

Сравним потенциальную энергию связи и кинетическую энергию для каждой из 5 систем.

In []: 
%%bash
for i in {'be','vr','nh','an','sd'}; do 
g_energy -f et_${i}.edr -o et_${i}_en.xvg -xvg none
done

Получили пять файлов в тремя столбиками время, потенциальная и кинетическая энергии

In []: 
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
In [83]: 
a = np.loadtxt("et_be_en.xvg")
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x=a[:,0]
y=a[:,1]
z=a[:,2]
print y, z
ax.scatter(x, y, z, alpha=0.5)

ax.set_xlabel('time, ps')
ax.set_ylabel('Potential, kJ/mol')
ax.set_zlabel('Kinetic, kJ/mol')

plt.show()

[  1.324735   4.551724  21.209444 ...,  27.335781  27.964699  26.9627  ] [  0.080536   5.043224   5.55969  ...,  24.721302  23.48753   24.736101]
In [85]: 
a = np.loadtxt("et_be_en.xvg")
x=a[:,0]
y=a[:,1]
print y
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('Potential, kJ/mol')

plt.show()

[  1.324735   4.551724  21.209444 ...,  27.335781  27.964699  26.9627  ]
In [86]: 
a = np.loadtxt("et_be_en.xvg")
x=a[:,0]
y=a[:,2]
print y
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('Kinetic, kJ/mol')

plt.show()

[  0.080536   5.043224   5.55969  ...,  24.721302  23.48753   24.736101]
In [87]: 
a = np.loadtxt("et_vr_en.xvg")
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x=a[:,0]
y=a[:,1]
z=a[:,2]
print y, z
ax.scatter(x, y, z, alpha=0.5)

ax.set_xlabel('time, ps')
ax.set_ylabel('Potential, kJ/mol')
ax.set_zlabel('Kinetic, kJ/mol')

plt.show()

[  1.324735  11.345387  16.668554 ...,  30.536285  21.279123  18.934366] [  0.080536  20.955555  18.634758 ...,  14.130837  11.77151   14.440428]
In [89]: 
a = np.loadtxt("et_be_en.xvg")
x=a[:,0]
y=a[:,2]
print y
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('Kinetic, kJ/mol')

plt.show()

[  0.080536   5.043224   5.55969  ...,  24.721302  23.48753   24.736101]
In [90]: 
a = np.loadtxt("et_nh_en.xvg")
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x=a[:,0]
y=a[:,1]
z=a[:,2]
print y, z
ax.scatter(x, y, z, alpha=0.5)

ax.set_xlabel('time, ps')
ax.set_ylabel('Potential, kJ/mol')
ax.set_zlabel('Kinetic, kJ/mol')

plt.show()

[  1.324735  54.647049   3.639566 ...,   6.11616    9.643476  29.092331] [  8.05360000e-02   1.00327408e+02   1.65455700e+00 ...,   3.30284880e+01
   1.21131490e+01   8.76123200e+01]
In [91]: 
a = np.loadtxt("et_an_en.xvg")
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x=a[:,0]
y=a[:,1]
z=a[:,2]
print y, z
ax.scatter(x, y, z, alpha=0.5)

ax.set_xlabel('time, ps')
ax.set_ylabel('Potential, kJ/mol')
ax.set_zlabel('Kinetic, kJ/mol')

plt.show()

[ 1.324735  0.528498  0.838018 ...,  0.577402  0.828243  1.026768] [ 0.080536  0.818427  0.536271 ...,  0.813213  0.552049  0.380507]
In [92]: 
a = np.loadtxt("et_an_en.xvg")
x=a[:,0]
y=a[:,1]
print y
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('Potential, kJ/mol')

plt.show()

[ 1.324735  0.528498  0.838018 ...,  0.577402  0.828243  1.026768]
In [93]: 
a = np.loadtxt("et_an_en.xvg")
x=a[:,0]
y=a[:,2]
print y
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('Kinetic, kJ/mol')

plt.show()

[ 0.080536  0.818427  0.536271 ...,  0.813213  0.552049  0.380507]
In [94]: 
a = np.loadtxt("et_sd_en.xvg")
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x=a[:,0]
y=a[:,1]
z=a[:,2]
print y, z
ax.scatter(x, y, z, alpha=0.5)

ax.set_xlabel('time, ps')
ax.set_ylabel('Potential, kJ/mol')
ax.set_zlabel('Kinetic, kJ/mol')

plt.show()

[  1.324735  11.847941  15.090025 ...,  21.228107  22.516361  22.121052] [  0.147808   7.453035  15.644512 ...,  24.095325  34.308815  31.656128]

Рассмотрим распределение длинны связи С-С за время моделирования

Создадим индекс файл b.ndx с одной связью в текстовом редакторе и запустим утилиту по анализу связей g_bond:

In []: 
%%bash
for i in {'be','vr','nh','an','sd'}; do 
g_bond -f et_${i}.trr -s et_${i}.tpr -o bond_${i}.xvg -n b.ndx -xvg none
done
In [127]: 
a = np.loadtxt("distance.xvg")#sp


plt.grid(True)


plt.subplot(211)
x=a[:,1]


n, bins, patches = plt.hist(x, 50, facecolor='m', alpha=0.75)

plt.xlabel('Distance')

plt.subplot(212)
x=a[:,0]
y=a[:,1]
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('dist')


plt.show()
In [126]: 
a = np.loadtxt("#distance.xvg.1#")#be


plt.grid(True)


plt.subplot(211)
x=a[:,1]


n, bins, patches = plt.hist(x, 50, facecolor='m', alpha=0.75)

plt.xlabel('Distance')

plt.subplot(212)
x=a[:,0]
y=a[:,1]
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('dist')


plt.show()
In [125]: 
a = np.loadtxt("#distance.xvg.2#")#vr


plt.grid(True)


plt.subplot(211)
x=a[:,1]


n, bins, patches = plt.hist(x, 50, facecolor='m', alpha=0.75)

plt.xlabel('Distance')

plt.subplot(212)
x=a[:,0]
y=a[:,1]
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('dist')


plt.show()
In [124]: 
a = np.loadtxt("#distance.xvg.3#")#nh


plt.grid(True)


plt.subplot(211)
x=a[:,1]


n, bins, patches = plt.hist(x, 50, facecolor='m', alpha=0.75)

plt.xlabel('Distance')

plt.subplot(212)
x=a[:,0]
y=a[:,1]
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('dist')


plt.show()
In [123]: 
a = np.loadtxt("#distance.xvg.4#")#an


plt.grid(True)


plt.subplot(211)
x=a[:,1]


n, bins, patches = plt.hist(x, 50, facecolor='m', alpha=0.75)

plt.xlabel('Distance')

plt.subplot(212)
x=a[:,0]
y=a[:,1]
plt.scatter(x, y)

plt.xlabel('time, ps')
plt.ylabel('dist')


plt.show()