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

Изучение работы методов контроля температуры в молекулярной динамике на примере моделирования молекулы этана с помощью алгоритмов метода Берендсена для контроля температуры, метода "Velocity rescale" для контроля температуры, метода Нуза-Хувера для контроля температуры, метода Андерсена для контроля температуры, метода стохастической молекулярной динамики и программы GROMACS.

In [ ]:
Дан файл с координатами атомов этана (http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro).
Дан файл с топологией этана, в нём надо заменить тип атомов, дописать связи и уголы второго атома C аналогично первому.
    Изначальный файл:
#include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp"

[ moleculetype ]
; Name            nrexcl
et            3

[ atoms ]
;   nr  type  resnr  residue  atom   cgnr     charge       mass
    1   CX      1    ETH      C1      1    -0.189      12.01
    2   CX      1    ETH      C2      2    -0.155      12.01
    3   HX      1    ETH      H1      3     0.0059       1.008
    4   HX      1    ETH      H2      4     0.0059       1.008
    5   HX      1    ETH      H3      5     0.0059       1.008
    6   HX      1    ETH      H4      6     0.0056       1.008
    7   HX      1    ETH      H5      7     0.0056       1.008
    8   HX      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  
........... 
надо дописать

[ 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
........... 
надо дописать 

[ dihedrals ]
;  ai    aj    ak    al funct  
    3    1     2     6      1  
    3    1     2     7      1 
    3    1     2     8      1  
    4    1     2     6      1  
 ........... 
надо дописать 
 
[ 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
    
    
    
   Исправленный файл:
#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
    
[ 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     8     1  
    1     2     7     1  
    1     2     6     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
sad
[ molecules ]
;Name count
 et    1

Далее на kodomo.¶

In [ ]:
Grompp - препроцессор gromacs, читает файл топологии молекул(ы), координат и .mdp-файл с алгоритмом контроля температуры, которых у нас пять:
be.mdp - метод Берендсена для контроля температуры.
vr.mdp - метод "Velocity rescale" для контроля температуры.
nh.mdp - метод Нуза-Хувера для контроля температуры.
an.mdp - метод Андерсена для контроля температуры.
sd.mdp - метод стохастической молекулярной динамики.

-f (grompp input file with Molecular Dinamics Parameters)
-c (coordinate file)
-p (topology file)
-o (output)  

Команда:
  gmx grompp -f be.mdp -c etane.gro -p et1.top -o etane_be.tpr

Ответ:

Command line:                
  gmx grompp -f be.mdp -c etane.gro -p et1.top -o etane_be.tpr

Ignoring obsolete mdp entry 'title'
Ignoring obsolete mdp entry 'cpp'
Ignoring obsolete mdp entry 'ns_type'
Replacing old mdp entry 'unconstrained-start' by 'continuation'

NOTE 1 [file be.mdp]:
  nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting
  nstcomm to nstcalcenergy


NOTE 2 [file be.mdp]:
  The Berendsen thermostat does not generate the correct kinetic energy
  distribution. You might want to consider using the V-rescale thermostat.

Setting the LD random seed to 813682169
Generated 330891 of the 330891 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 330891 of the 330891 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'et'

NOTE 3 [file et1.top, line 77]:
  System has non-zero total charge: -0.305500
  Total charge should normally be an integer. See
  http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
  for discussion on how close it should be to an integer.



Analysing residue names:
There are:     1      Other residues
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Number of degrees of freedom in T-Coupling group System is 21.00
Determining Verlet buffer for a tolerance of 0.005 kJ/mol/ps at 300 K
Calculated rlist for 1x1 atom pair-list as 0.702 nm, buffer size 0.002 nm
Set rlist, assuming 4x4 atom pair-list, to 0.700 nm, buffer size 0.000 nm
Note that mdrun will redetermine rlist based on the actual pair-list setup

NOTE 4 [file be.mdp]:
  You are using a plain Coulomb cut-off, which might produce artifacts.
  You might want to consider using PME electrostatics.


This run will generate roughly 10 Mb of data

There were 4 notes

GROMACS reminds you: "Can someone please tell Icarus that he's not the only one falling from the sky?" (Urban Dance Squad)
In [ ]:
Затем используем mdrun - главный вычислительный двигатель gromacs'a.

-deffnm - имя для выходных файлов
-nt - сколько потоков процессора использовать (?)

Пишет логи в один файл (-g). Пишет файл траекторий (-o), содержащий координаты, скорости и дополнительные силы. Пишет структурный файл (-c), содержащий координаты и скорости последнего этапа. Пишет файл энергий (-e), содержащий энергии, температуру, давление и проч.

Команда: 
   gmx mdrun -deffnm etane_be -v -nt 1

Ответ:

Command line:
  gmx mdrun -deffnm etane_be -v -nt 1

Compiled SIMD: SSE4.1, but for this host/run AVX2_256 might be better (see
log).
The current CPU can measure timings more accurately than the code in
gmx mdrun was configured to use. This might affect your simulation
speed as accurate timings are needed for load-balancing.
Please consider rebuilding gmx mdrun with the GMX_USE_RDTSCP=ON CMake option.
Reading file etane_be.tpr, VERSION 2020.6-Debian-2020.6-2 (single precision)

NOTE: Parallelization is limited by the small number of atoms,
      only starting 1 thread-MPI ranks.
      You can use the -nt and/or -ntmpi option to optimize the number of threads.

Changing nstlist from 10 to 50, rlist from 0.7 to 0.746


Using 1 MPI thread

Non-default thread affinity set, disabling internal thread affinity

Using 1 OpenMP thread

starting mdrun 'sad'
250000 steps,    250.0 ps.
step 249900, remaining wall clock time:     0 s
Writing final coordinates.
step 250000, remaining wall clock time:     0 s
               Core t (s)   Wall t (s)        (%)
       Time:        1.800        1.800      100.0
                 (ns/day)    (hour/ns)
Performance:    11998.437        0.002

GROMACS reminds you: "What If None Of Your Dreams Come True ?" (E. Costello)
Я бы предпочла, чтобы он этого не делал.
In [ ]:
Из файла с посчитанными энергиями найдём нужные командой energy.

-f - файл с энергиями 
-o - выходной файл
-xvg - график

Команда: 
    gmx energy -f etane_be.edr -o etane_be_energy_pot.xvg -xvg none

Ответ:
    
Opened etane_be.edr as single precision energy file

Select the terms you want from the following list by
selecting either (part of) the name or the number or a combination.
End your selection with an empty line or a zero.
-------------------------------------------------------------------
  1  Bond             2  Angle            3  Ryckaert-Bell.   4  LJ-14
  5  Coulomb-14       6  LJ-(SR)          7  Coulomb-(SR)     8  Potential
  9  Kinetic-En.     10  Total-Energy    11  Conserved-En.   12  Temperature
 13  Pressure        14  Vir-XX          15  Vir-XY          16  Vir-XZ
 17  Vir-YX          18  Vir-YY          19  Vir-YZ          20  Vir-ZX
 21  Vir-ZY          22  Vir-ZZ          23  Pres-XX         24  Pres-XY
 25  Pres-XZ         26  Pres-YX         27  Pres-YY         28  Pres-YZ
 29  Pres-ZX         30  Pres-ZY         31  Pres-ZZ         32  #Surf*SurfTen
 33  T-System        34  Lamb-System

8

Last energy frame read 2500 time  250.000

Statistics over 250001 steps [ 0.0000 through 250.0000 ps ], 1 data sets
All statistics are over 2501 points (frames)

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Potential                   17.0246      0.086    1.62639   0.542946  (kJ/mol)

GROMACS reminds you: "The absence of real intelligence doesn't prove you're using AI" (Magnus Lundborg)

    
Команда: 
    gmx energy -f etane_be.edr -o etane_be_energy_kin.xvg -xvg none
    
Ответ:

9


Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Kinetic En.                 26.1814     0.0068    1.60127 0.00367745  (kJ/mol)

GROMACS reminds you: "In ancient times they had no statistics so they had to fall back on lies." (Stephen Leacock)
    
Команда: 
    gmx energy -f etane_be.edr -o etane_be_energy_total.xvg -xvg none
    
Ответ:
    
10


Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy                43.2061      0.087    1.47225   0.546623  (kJ/mol)

GROMACS reminds you: "The Universe is Somewhere In Here" (J.G.E.M. Fraaije)
In [ ]:
Затем запустим trjconv, которая конвертирует файлы траекторий. Нам нужны PDB-файлы для перехода к виртуальной репрезентации моделирвоания.

-f - файл траекторий
-s - дополнительный файл, например, выход первой подпрограммы с топологией и параметрами симуляции
-o - выходной файл


Команда: 
    gmx trjconv -f etane_be.trr -s etane_be.tpr -o etane_be.pdb
    
Ответ (без лишних частей):

    Command line:
  gmx trjconv -f etane_be.trr -s etane_be.tpr -o etane_be.pdb

Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write pdb: Protein data bank file
Reading file etane_be.tpr, VERSION 2020.6-Debian-2020.6-2 (single precision)
Reading file etane_be.tpr, VERSION 2020.6-Debian-2020.6-2 (single precision)
Select group for output
Group     0 (         System) has     8 elements
Group     1 (          Other) has     8 elements
Group     2 (            ETH) has     8 elements
Select a group: 0
Selected 0: 'System'
trr version: GMX_trn_file (single precision)
 ->  frame    250 time  250.000        ->  frame    200 time  200.000


GROMACS reminds you: "Where all think alike, no one thinks very much." (Walter Lippmann)

Command line:
  gmx energy -f etane_be.edr -o etane_be_bond.xvg -xvg none

Opened etane_be.edr as double precision energy file

Select the terms you want from the following list by
selecting either (part of) the name or the number or a combination.
End your selection with an empty line or a zero.
-------------------------------------------------------------------
  1  Bond             2  Angle            3  Ryckaert-Bell.   4  LJ-14
  5  Coulomb-14       6  LJ-(SR)          7  Coulomb-(SR)     8  Potential
  9  Kinetic-En.     10  Total-Energy    11  Conserved-En.   12  Temperature
 13  Pressure        14  Vir-XX          15  Vir-XY          16  Vir-XZ
 17  Vir-YX          18  Vir-YY          19  Vir-YZ          20  Vir-ZX
 21  Vir-ZY          22  Vir-ZZ          23  Pres-XX         24  Pres-XY
 25  Pres-XZ         26  Pres-YX         27  Pres-YY         28  Pres-YZ
 29  Pres-ZX         30  Pres-ZY         31  Pres-ZZ         32  #Surf*SurfTen
 33  T-System        34  Lamb-System

1

Last energy frame read 2500 time  250.000

Statistics over 250001 steps [ 0.0000 through 250.0000 ps ], 1 data sets
All statistics are over 2501 points (frames)

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Bond                        2.23315       0.88    2.84702   -5.76501  (kJ/mol)

GROMACS reminds you: "Contemplating answers that could break my bonds." (Peter Hammill)

Визуальный анализ¶

In [1]:
import numpy as np
from IPython.display import Image

import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd, stored
In [9]:
import os
import imageio
import time
from IPython.display import Image
In [3]:
thermo_types = ['be', 'vr', 'nh', 'sd']
for thermo in thermo_types:
    # aesthetics
    cmd.reinitialize()
    cmd.set('sphere_scale', 0.2)
    cmd.set('stick_radius', 0.1)
    cmd.set('ray_trace_mode', 3)
    cmd.set('antialias', 2)
    cmd.set('sphere_scale', 0.2)
    cmd.set('stick_radius', 0.1)
    # load data
    cmd.load(f'D:\pain_pr7\etane_{thermo}.pdb')
    cmd.show('spheres')
    cmd.orient()
    # save images
    os.mkdir(thermo)
    cmd.mpng(os.path.join(thermo, ''), mode=2, width=640, height=480)
    while not os.path.exists(os.path.join(thermo, '0251.png')):
        time.sleep(0.1)
    # create gif    
    with imageio.get_writer(f'{thermo}.gif', mode='I', fps=30) as writer:
        for i in range(0, 251):
            image = imageio.imread(os.path.join(thermo, f'{i + 1:04}.png'))
            writer.append_data(image)

Моделирование с методом Берендсена для контроля температуры:

In [10]:
Image(filename='be.gif')
Out[10]:
<IPython.core.display.Image object>

Из-за потери колебательных степеней свобод атомов мы не наблюдаем заслоненной конформации, CH3-группы равномерно вращаются, их колебания перешли в поступательное движение туда-сюда.

Моделирование с методом "Velocity rescale" для контроля температуры:

In [11]:
Image(filename='vr.gif')
Out[11]:
<IPython.core.display.Image object>

Моделирование с методом Нуза-Хувера для контроля температуры:

In [12]:
Image(filename='nh.gif')
Out[12]:
<IPython.core.display.Image object>

Моделирование с методом стохастической молекулярной динамики:

In [14]:
Image(filename='sd.gif')
Out[14]:
<IPython.core.display.Image object>
In [ ]:
Моделирование с методом Андерсена для контроля температуры:
Для Андерсона не считает.
    
Command line:
  gmx grompp -f an.mdp -c etane.gro -p et1.top -o etane_an.tpr

Ignoring obsolete mdp entry 'title'
Ignoring obsolete mdp entry 'cpp'
Ignoring obsolete mdp entry 'ns_type'
Replacing old mdp entry 'unconstrained-start' by 'continuation'

NOTE 1 [file an.mdp]:
  nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting
  nstcomm to nstcalcenergy


ERROR 1 [file an.mdp]:
  Andersen temperature control not supported for integrator md.


NOTE 2 [file an.mdp]:
  Center of mass removal not necessary for Andersen.  All velocities of
  coupled groups are rerandomized periodically, so flying ice cube errors
  will not occur.


ERROR 2 [file an.mdp]:
  nstcomm must be 1, not 100 for Andersen, as velocities of atoms in
  coupled groups are randomized every time step

Setting the LD random seed to -19542161
Generated 330891 of the 330891 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 330891 of the 330891 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'et'

NOTE 3 [file et1.top, line 77]:
  System has non-zero total charge: -0.305500
  Total charge should normally be an integer. See
  http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
  for discussion on how close it should be to an integer.




There were 3 notes

-------------------------------------------------------
Program:     gmx grompp, version 2020.6-Debian-2020.6-2
Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 1928)

Fatal error:
There were 2 errors in input file(s)

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------

Графики изменения энергий

Объединяем файлы энергий в одном файле, строим графики.

In [ ]:
Метод Берендсена
In [22]:
import matplotlib.pyplot as plt
time, pot, kin, tot = [], [], [], []

for line in open('etane_be_energy.txt', 'r'):
    values = [float(s) for s in line.split('\t')]
    time.append(values[0])
    pot.append(values[1])
    kin.append(values[2])    
    tot.append(values[3])

plt.plot(time, pot, label='potential')
plt.plot(time, kin, label='kinetic')
plt.plot(time, tot, label='total')
plt.legend(title='Energies', loc='lower right')

plt.xlabel('Time, ps')
plt.ylabel('Energy, kJ/mol')

plt.show()
In [ ]:
Метод "Velocity rescale"
In [27]:
import matplotlib.pyplot as plt
time, pot, kin, tot = [], [], [], []

for line in open('etane_vr_energy.txt', 'r'):
    values = [float(s) for s in line.split('\t')]
    time.append(values[0])
    pot.append(values[1])
    kin.append(values[2])    
    tot.append(values[3])

plt.plot(time, pot, label='potential')
plt.plot(time, kin, label='kinetic')
plt.plot(time, tot, label='total')
plt.legend(title='Energies', loc='lower right')

plt.xlabel('Time, ps')
plt.ylabel('Energy, kJ/mol')

plt.show()
In [ ]:
Метод Нуза-Хувера
In [25]:
import matplotlib.pyplot as plt
time, pot, kin, tot = [], [], [], []

for line in open('etane_nh_energy.txt', 'r'):
    values = [float(s) for s in line.split('\t')]
    time.append(values[0])
    pot.append(values[1])
    kin.append(values[2])    
    tot.append(values[3])

plt.plot(time, pot, label='potential')
plt.plot(time, kin, label='kinetic')
plt.plot(time, tot, label='total')
plt.legend(title='Energies', loc='lower right')

plt.xlabel('Time, ps')
plt.ylabel('Energy, kJ/mol')

plt.show()

Метод стохастической молекулярной динамики

In [28]:
import matplotlib.pyplot as plt
time, pot, kin, tot = [], [], [], []

for line in open('etane_sd_energy.txt', 'r'):
    values = [float(s) for s in line.split('\t')]
    time.append(values[0])
    pot.append(values[1])
    kin.append(values[2])    
    tot.append(values[3])

plt.plot(time, pot, label='potential')
plt.plot(time, kin, label='kinetic')
plt.plot(time, tot, label='total')
plt.legend(title='Energies', loc='lower right')

plt.xlabel('Time, ps')
plt.ylabel('Energy, kJ/mol')

plt.show()

Лучше всего выглядят velocity rescale и метод стохастической молекулярной динамики, а гифка лучшая вышла у velocity rescale.