Подготовка системы: et.gro, et.top, mdp-файлы
In [1]:
import sys
print(sys.executable)
!which gmx
!gmx --version
/Users/romashka/Documents/jupyter/envs/jlab/bin/python
/Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
:-) GROMACS - gmx, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx --version
GROMACS version: 2026.0-conda_forge
Precision: mixed
Memory model: 64 bit
MPI library: thread_mpi
MPI version: built in
OpenMP support: enabled (GMX_OPENMP_MAX_THREADS = 128)
GPU support: disabled
SIMD instructions: ARM_NEON_ASIMD
CPU FFT library: fftw-3.3.10
GPU FFT library: none
Multi-GPU FFT: none
TNG support: enabled
Hwloc support: disabled
Tracing support: disabled
Colvars support: enabled (version 2025-10-13)
CP2K support: disabled
Torch support: disabled
Plumed support: enabled
C compiler: /Users/runner/miniforge3/conda-bld/gromacs_1770909788555/_build_env/bin/arm64-apple-darwin20.0.0-clang Clang 19.1.7
C compiler flags: -Wno-missing-field-initializers -O3 -DNDEBUG
C++ compiler: /Users/runner/miniforge3/conda-bld/gromacs_1770909788555/_build_env/bin/arm64-apple-darwin20.0.0-clang++ Clang 19.1.7
C++ compiler flags: -Wno-reserved-identifier -Wno-missing-field-initializers -Weverything -Wno-c++98-compat -Wno-c++98-compat-pedantic -Wno-source-uses-openmp -Wno-c++17-extensions -Wno-documentation-unknown-command -Wno-covered-switch-default -Wno-switch-enum -Wno-switch-default -Wno-extra-semi-stmt -Wno-weak-vtables -Wno-shadow -Wno-padded -Wno-reserved-id-macro -Wno-double-promotion -Wno-exit-time-destructors -Wno-global-constructors -Wno-documentation -Wno-format-nonliteral -Wno-used-but-marked-unused -Wno-float-equal -Wno-conditional-uninitialized -Wno-conversion -Wno-disabled-macro-expansion -Wno-unused-macros -Wno-unsafe-buffer-usage -Wno-cast-function-type-strict SHELL:-fopenmp=libomp -O3 -DNDEBUG
BLAS library: External - detected on the system
LAPACK library: External - detected on the system
In [ ]:
In [2]:
import os
workdir = "/Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats"
os.makedirs(workdir, exist_ok=True)
os.chdir(workdir)
print("Рабочая папка:")
print(os.getcwd())
Рабочая папка: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
In [3]:
!curl -L -o et.gro http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
print("Файлы в папке:")
!ls -lh
print("Содержимое et.gro:")
!cat et.gro
% Total % Received % Xferd Average Speed Time Time Time Current
Dload Upload Total Spent Left Speed
100 384 100 384 0 0 1235 0 --:--:-- --:--:-- --:--:-- 1238
100 399 100 399 0 0 619 0 --:--:-- --:--:-- --:--:-- 619
Файлы в папке:
total 24
-rw-r--r-- 1 romashka staff 399B 27 мая 19:35 et.gro
-rw-r--r-- 1 romashka staff 4,1K 27 мая 19:35 task7_ethane_thermostats.ipynb
Содержимое et.gro:
etane
8
1ETH C1 1 0.577 0.217 0.574
1ETH C2 2 0.680 0.252 0.467
1ETH H1 3 0.478 0.241 0.538
1ETH H2 4 0.597 0.274 0.664
1ETH H3 5 0.583 0.111 0.597
1ETH H4 6 0.676 0.358 0.445
1ETH H5 7 0.660 0.195 0.377
1ETH H6 8 0.780 0.228 0.504
1.50000 1.50000 1.50000
In [4]:
print("Тип углерода для CH3:")
!grep -n "opls_135" /Users/romashka/Documents/jupyter/envs/jlab/share/gromacs/top/oplsaa.ff/atomtypes.atp
print("\nТип водорода алкана:")
!grep -n "opls_140" /Users/romashka/Documents/jupyter/envs/jlab/share/gromacs/top/oplsaa.ff/atomtypes.atp
Тип углерода для CH3: 8:; Explicit all-atom parameters start with opls_135. 148: opls_135 12.01100 ; alkane CH3 Тип водорода алкана: 153: opls_140 1.00800 ; alkane H.
In [5]:
!grep -n "HC.*CT.*CT.*HC" /Users/romashka/Documents/jupyter/envs/jlab/share/gromacs/top/oplsaa.ff/ffbonded.itp | head
1119: CT CO HC 1 110.700 313.800 ; -idem- : CT-CT-HC- wd 6/95 Glucose 1123: CO CT HC 1 110.700 313.800 ; -idem- : CT-CT-HC- wd 6/95 Glucose 1824: HC CT CT HC 3 0.62760 1.88280 0.00000 -2.51040 0.00000 0.00000 ; hydrocarbon *new* 11/99 1855: HC CT CT_2 HC 3 0.62760 1.88280 0.00000 -2.51040 0.00000 0.00000 ; hydrocarbon all-atom 1861: HC CT_3 CT HC 3 0.62760 1.88280 0.00000 -2.51040 0.00000 0.00000 ;
In [6]:
top_text = """
#include "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.180 12.011
2 opls_135 1 ETH C2 2 -0.180 12.011
3 opls_140 1 ETH H1 3 0.060 1.008
4 opls_140 1 ETH H2 4 0.060 1.008
5 opls_140 1 ETH H3 5 0.060 1.008
6 opls_140 1 ETH H4 6 0.060 1.008
7 opls_140 1 ETH H5 7 0.060 1.008
8 opls_140 1 ETH H6 8 0.060 1.008
[ bonds ]
; ai aj funct
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
; around C1
3 1 4 1
3 1 5 1
4 1 5 1
2 1 3 1
2 1 4 1
2 1 5 1
; around C2
6 2 7 1
6 2 8 1
7 2 8 1
1 2 6 1
1 2 7 1
1 2 8 1
[ dihedrals ]
; ai aj ak al funct
; OPLS torsions use funct = 3
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 ]
; ai aj funct
3 6 1
3 7 1
3 8 1
4 6 1
4 7 1
4 8 1
5 6 1
5 7 1
5 8 1
[ system ]
one ethane molecule
[ molecules ]
; Name count
et 1
"""
with open("et.top", "w") as f:
f.write(top_text.strip() + "\n")
print("Файл et.top создан")
Файл et.top создан
In [7]:
!ls -lh
!head -60 et.top
!ls et.gro et.top
total 32 -rw-r--r-- 1 romashka staff 399B 27 мая 19:35 et.gro -rw-r--r-- 1 romashka staff 1,5K 27 мая 19:37 et.top -rw-r--r-- 1 romashka staff 7,5K 27 мая 19:37 task7_ethane_thermostats.ipynb #include "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.180 12.011 2 opls_135 1 ETH C2 2 -0.180 12.011 3 opls_140 1 ETH H1 3 0.060 1.008 4 opls_140 1 ETH H2 4 0.060 1.008 5 opls_140 1 ETH H3 5 0.060 1.008 6 opls_140 1 ETH H4 6 0.060 1.008 7 opls_140 1 ETH H5 7 0.060 1.008 8 opls_140 1 ETH H6 8 0.060 1.008 [ bonds ] ; ai aj funct 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 ; around C1 3 1 4 1 3 1 5 1 4 1 5 1 2 1 3 1 2 1 4 1 2 1 5 1 ; around C2 6 2 7 1 6 2 8 1 7 2 8 1 1 2 6 1 1 2 7 1 1 2 8 1 [ dihedrals ] ; ai aj ak al funct ; OPLS torsions use funct = 3 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 ] ; ai aj funct et.gro et.top
In [9]:
mdp_common = """
; common molecular dynamics settings for ethane
integrator = md-vv
dt = 0.001
nsteps = 50000
; output
nstxout = 10
nstvout = 10
nstfout = 0
nstlog = 100
nstenergy = 10
nstcalcenergy = 10
; no constraints, because we want to analyze C-C bond length
constraints = none
; nonbonded settings
cutoff-scheme = Verlet
nstlist = 10
rlist = 1.0
coulombtype = Cut-off
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
; one molecule in vacuum
pbc = no
; center of mass motion
comm-mode = Linear
nstcomm = 100
; initial velocities
gen-vel = yes
gen-temp = 320
gen-seed = 7777
; temperature group
tc-grps = System
ref-t = 320
"""
mdp_files = {
"be.mdp": mdp_common + """
; Berendsen thermostat
tcoupl = Berendsen
tau-t = 0.1
""",
"vr.mdp": mdp_common + """
; Velocity rescale thermostat
tcoupl = V-rescale
tau-t = 0.1
""",
"nh.mdp": mdp_common + """
; Nose-Hoover thermostat
tcoupl = Nose-Hoover
tau-t = 0.5
""",
"an.mdp": mdp_common + """
; Andersen thermostat
tcoupl = Andersen
tau-t = 0.1
""",
"sd.mdp": """
; stochastic dynamics for ethane
integrator = sd
dt = 0.001
nsteps = 50000
; output
nstxout = 10
nstvout = 10
nstfout = 0
nstlog = 100
nstenergy = 10
nstcalcenergy = 10
constraints = none
cutoff-scheme = Verlet
nstlist = 10
rlist = 1.0
coulombtype = Cut-off
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
pbc = no
comm-mode = Linear
nstcomm = 100
; stochastic dynamics temperature settings
tc-grps = System
tau-t = 0.1
ref-t = 320
gen-vel = yes
gen-temp = 320
gen-seed = 7777
ld-seed = 7777
"""
}
for name, text in mdp_files.items():
with open(name, "w") as f:
f.write(text.strip() + "\n")
print("Созданы файлы:")
!ls -lh *.mdp
Созданы файлы: -rw-r--r-- 1 romashka staff 1,1K 27 мая 19:40 an.mdp -rw-r--r-- 1 romashka staff 1,1K 27 мая 19:40 be.mdp -rw-r--r-- 1 romashka staff 1,1K 27 мая 19:40 nh.mdp -rw-r--r-- 1 romashka staff 918B 27 мая 19:40 sd.mdp -rw-r--r-- 1 romashka staff 1,1K 27 мая 19:40 vr.mdp
In [10]:
!cat be.mdp
; common molecular dynamics settings for ethane integrator = md-vv dt = 0.001 nsteps = 50000 ; output nstxout = 10 nstvout = 10 nstfout = 0 nstlog = 100 nstenergy = 10 nstcalcenergy = 10 ; no constraints, because we want to analyze C-C bond length constraints = none ; nonbonded settings cutoff-scheme = Verlet nstlist = 10 rlist = 1.0 coulombtype = Cut-off rcoulomb = 1.0 vdwtype = Cut-off rvdw = 1.0 ; one molecule in vacuum pbc = no ; center of mass motion comm-mode = Linear nstcomm = 100 ; initial velocities gen-vel = yes gen-temp = 320 gen-seed = 7777 ; temperature group tc-grps = System ref-t = 320 ; Berendsen thermostat tcoupl = Berendsen tau-t = 0.1
Создаем коробку для этана
In [13]:
!gmx editconf -f et.gro -o et_box.gro -c -d 1.5 -bt cubic
:-) GROMACS - gmx editconf, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx editconf -f et.gro -o et_box.gro -c -d 1.5 -bt cubic
Note that major changes are planned in future for editconf, to improve usability and utility.
Read 8 atoms
Volume: 3.375 nm^3, corresponds to roughly 1500 electrons
No velocities found
system size : 0.302 0.247 0.287 (nm)
diameter : 0.305 (nm)
center : 0.629 0.235 0.521 (nm)
box vectors : 1.500 1.500 1.500 (nm)
box angles : 90.00 90.00 90.00 (degrees)
box volume : 3.38 (nm^3)
shift : 1.023 1.418 1.132 (nm)
new center : 1.652 1.652 1.652 (nm)
new box vectors : 3.305 3.305 3.305 (nm)
new box angles : 90.00 90.00 90.00 (degrees)
new box volume : 36.09 (nm^3)
GROMACS reminds you: "Do you know what cations don't like? Dog-ions. Do you know what they like? Pie." (Tom Cheatham)
In [14]:
!ls -lh et.gro et_box.gro
!cat et_box.gro
-rw-r--r-- 1 romashka staff 403B 27 мая 19:44 et_box.gro
-rw-r--r-- 1 romashka staff 399B 27 мая 19:35 et.gro
etane
8
1ETH C1 1 1.600 1.635 1.706
1ETH C2 2 1.703 1.670 1.599
1ETH H1 3 1.501 1.659 1.670
1ETH H2 4 1.620 1.692 1.796
1ETH H3 5 1.606 1.529 1.729
1ETH H4 6 1.699 1.776 1.577
1ETH H5 7 1.683 1.613 1.509
1ETH H6 8 1.803 1.646 1.636
3.30457 3.30457 3.30457
Переписываем .mdp
In [15]:
for fname in ["be.mdp", "vr.mdp", "nh.mdp", "an.mdp", "sd.mdp"]:
with open(fname) as f:
text = f.read()
text = text.replace("pbc = no", "pbc = xyz")
if "verlet-buffer-tolerance" not in text:
text = text.replace(
"rvdw = 1.0",
"rvdw = 1.0\nverlet-buffer-tolerance = -1"
)
with open(fname, "w") as f:
f.write(text)
print("Исправила pbc и verlet-buffer-tolerance во всех mdp-файлах")
!grep -n "pbc\|verlet-buffer" *.mdp
Исправила pbc и verlet-buffer-tolerance во всех mdp-файлах an.mdp:26:verlet-buffer-tolerance = -1 an.mdp:29:pbc = xyz be.mdp:26:verlet-buffer-tolerance = -1 be.mdp:29:pbc = xyz mdout.mdp:89:pbc = no mdout.mdp:93:verlet-buffer-tolerance = 0.005 mdout.mdp:95:verlet-buffer-pressure-tolerance = 0.5 nh.mdp:26:verlet-buffer-tolerance = -1 nh.mdp:29:pbc = xyz sd.mdp:24:verlet-buffer-tolerance = -1 sd.mdp:26:pbc = xyz vr.mdp:26:verlet-buffer-tolerance = -1 vr.mdp:29:pbc = xyz
grompp для Берендсена
In [17]:
!gmx grompp -f be.mdp -c et_box.gro -p et.top -o et_be.tpr -maxwarn 1
:-) GROMACS - gmx grompp, 2026.0-conda_forge (-: Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx Data prefix: /Users/romashka/Documents/jupyter/envs/jlab Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats Command line: gmx grompp -f be.mdp -c et_box.gro -p et.top -o et_be.tpr -maxwarn 1 NOTE 1 [file be.mdp]: With Verlet lists and no buffer tolerance, the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation. NOTE 2 [file be.mdp]: rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation. NOTE 3 [file be.mdp]: verlet-buffer-pressure-tolerance is ignored when verlet-buffer-tolerance < 0 WARNING 1 [file be.mdp]: The Berendsen thermostat does not generate the correct kinetic energy distribution, and should not be used for new production simulations (in our opinion). We would recommend the V-rescale thermostat. Setting the LD random seed to -1897038373 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' Velocities were taken from a Maxwell distribution at 320 K 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 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 8 Mb of data There were 4 NOTEs There was 1 WARNING GROMACS reminds you: "Inventions have long since reached their limit, and I see no hope for further development." (Julius Sextus Frontinus, 1st century A.D.)
In [18]:
!ls -lh et_be.tpr
-rw-r--r-- 1 romashka staff 5,4K 27 мая 19:47 et_be.tpr
запускаем первый mdrun
In [19]:
!gmx mdrun -deffnm et_be -v -nt 1
:-) GROMACS - gmx mdrun, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx mdrun -deffnm et_be -v -nt 1
Reading file et_be.tpr, VERSION 2026.0-conda_forge (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.
Can not increase nstlist because verlet-buffer-tolerance is not set or used
Using 1 MPI thread
Using 1 OpenMP thread
starting mdrun 'one ethane molecule'
50000 steps, 50.0 ps.
step 50000, remaining wall clock time: 0 s
Writing final coordinates.
Core t (s) Wall t (s) (%)
Time: 0.138 0.138 100.0
(ns/day) (hour/ns) (ms/step) (Matom*steps/s)
Performance: 31286.161 0.001 0.003 2.897
GROMACS reminds you: "For those who want some proof that physicists are human, the proof is in the idiocy of all the different units which they use for measuring energy." (Richard Feynman)
In [20]:
!ls -lh et_be.*
-rw-r--r-- 1 romashka staff 3,1K 27 мая 19:48 et_be.cpt -rw-r--r-- 1 romashka staff 1,0M 27 мая 19:48 et_be.edr -rw-r--r-- 1 romashka staff 609B 27 мая 19:48 et_be.gro -rw-r--r-- 1 romashka staff 255K 27 мая 19:48 et_be.log -rw-r--r-- 1 romashka staff 5,4K 27 мая 19:47 et_be.tpr -rw-r--r-- 1 romashka staff 1,5M 27 мая 19:48 et_be.trr
Делаем .tpr для остальных методов
In [21]:
methods = ["vr", "nh", "an", "sd"]
for m in methods:
print("=" * 60)
print("Готовим:", m)
!gmx grompp -f {m}.mdp -c et_box.gro -p et.top -o et_{m}.tpr -maxwarn 1
============================================================
Готовим: vr
:-) GROMACS - gmx grompp, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx grompp -f vr.mdp -c et_box.gro -p et.top -o et_vr.tpr -maxwarn 1
NOTE 1 [file vr.mdp]:
With Verlet lists and no buffer tolerance, the optimal nstlist is >= 10,
with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect
on the accuracy of your simulation.
NOTE 2 [file vr.mdp]:
rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet
buffer. The cluster pair list does have a buffering effect, but choosing
a larger rlist might be necessary for good energy conservation.
NOTE 3 [file vr.mdp]:
verlet-buffer-pressure-tolerance is ignored when verlet-buffer-tolerance
< 0
Setting the LD random seed to -1140916269
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'
Velocities were taken from a Maxwell distribution at 320 K
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
NOTE 4 [file vr.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 8 Mb of data
There were 4 NOTEs
GROMACS reminds you: "Water is just water" (Berk Hess)
============================================================
Готовим: nh
:-) GROMACS - gmx grompp, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx grompp -f nh.mdp -c et_box.gro -p et.top -o et_nh.tpr -maxwarn 1
NOTE 1 [file nh.mdp]:
With Verlet lists and no buffer tolerance, the optimal nstlist is >= 10,
with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect
on the accuracy of your simulation.
NOTE 2 [file nh.mdp]:
rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet
buffer. The cluster pair list does have a buffering effect, but choosing
a larger rlist might be necessary for good energy conservation.
NOTE 3 [file nh.mdp]:
verlet-buffer-pressure-tolerance is ignored when verlet-buffer-tolerance
< 0
Setting the LD random seed to -1109133578
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'
Velocities were taken from a Maxwell distribution at 320 K
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
NOTE 4 [file nh.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 8 Mb of data
There were 4 NOTEs
GROMACS reminds you: "Water is just water" (Berk Hess)
============================================================
Готовим: an
:-) GROMACS - gmx grompp, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx grompp -f an.mdp -c et_box.gro -p et.top -o et_an.tpr -maxwarn 1
NOTE 1 [file an.mdp]:
With Verlet lists and no buffer tolerance, the optimal nstlist is >= 10,
with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect
on the accuracy of your simulation.
NOTE 2 [file an.mdp]:
rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet
buffer. The cluster pair list does have a buffering effect, but choosing
a larger rlist might be necessary for good energy conservation.
NOTE 3 [file an.mdp]:
verlet-buffer-pressure-tolerance is ignored when verlet-buffer-tolerance
< 0
NOTE 4 [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 1 [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 318485469
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'
Velocities were taken from a Maxwell distribution at 320 K
There were 4 NOTEs
-------------------------------------------------------
Program: gmx grompp, version 2026.0-conda_forge
Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 2413)
Fatal error:
There was 1 error in input file(s)
For more information and tips for troubleshooting, please check the GROMACS
website at https://manual.gromacs.org/current/user-guide/run-time-errors.html
-------------------------------------------------------
============================================================
Готовим: sd
:-) GROMACS - gmx grompp, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx grompp -f sd.mdp -c et_box.gro -p et.top -o et_sd.tpr -maxwarn 1
NOTE 1 [file sd.mdp]:
With Verlet lists and no buffer tolerance, the optimal nstlist is >= 10,
with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect
on the accuracy of your simulation.
NOTE 2 [file sd.mdp]:
rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet
buffer. The cluster pair list does have a buffering effect, but choosing
a larger rlist might be necessary for good energy conservation.
NOTE 3 [file sd.mdp]:
verlet-buffer-pressure-tolerance is ignored when verlet-buffer-tolerance
< 0
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'
Velocities were taken from a Maxwell distribution at 320 K
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
NOTE 4 [file sd.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 8 Mb of data
There were 4 NOTEs
GROMACS reminds you: "This Puke Stinks Like Beer" (LIVE)
In [22]:
!ls -lh *.tpr
-rw-r--r-- 1 romashka staff 5,4K 27 мая 19:47 et_be.tpr -rw-r--r-- 1 romashka staff 5,4K 27 мая 19:50 et_nh.tpr -rw-r--r-- 1 romashka staff 5,4K 27 мая 19:50 et_sd.tpr -rw-r--r-- 1 romashka staff 5,4K 27 мая 19:50 et_vr.tpr
In [25]:
with open("an.mdp") as f:
text = f.read()
text = text.replace("nstcomm = 100", "nstcomm = 1")
with open("an.mdp", "w") as f:
f.write(text)
print("Исправленный an.mdp:")
!grep -n "tcoupl\|nstcomm\|tau-t\|ref-t" an.mdp
Исправленный an.mdp: 33:nstcomm = 1 42:ref-t = 320 45:tcoupl = Andersen 46:tau-t = 0.1
In [26]:
!gmx grompp -f an.mdp -c et_box.gro -p et.top -o et_an.tpr -maxwarn 5
:-) GROMACS - gmx grompp, 2026.0-conda_forge (-: Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx Data prefix: /Users/romashka/Documents/jupyter/envs/jlab Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats Command line: gmx grompp -f an.mdp -c et_box.gro -p et.top -o et_an.tpr -maxwarn 5 NOTE 1 [file an.mdp]: With Verlet lists and no buffer tolerance, the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation. NOTE 2 [file an.mdp]: rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation. NOTE 3 [file an.mdp]: verlet-buffer-pressure-tolerance is ignored when verlet-buffer-tolerance < 0 NOTE 4 [file an.mdp]: nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, consider setting nstcomm equal to nstcalcenergy for less overhead NOTE 5 [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. Setting the LD random seed to -308332036 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' Velocities were taken from a Maxwell distribution at 320 K Analysing residue names: There are: 1 Other residues Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups... NOTE 6 [file an.mdp]: Andersen temperature control methods assume nsttcouple = 1; there is no need for larger nsttcouple > 1, since no global parameters are computed. nsttcouple has been reset to 1 Number of degrees of freedom in T-Coupling group System is 21.00 NOTE 7 [file an.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 8 Mb of data There were 7 NOTEs GROMACS reminds you: "Ludwig Boltzmann, who spent much of his life studying statistical mechanics, died in 1906, by his own hand. Paul Ehrenfest, carrying on the same work, died similarly in 1933. Now it is our turn to study statistical mechanics. Perhaps it will be wise to approach the subject cautiously." (David Goodstein)
In [27]:
!ls -lh et_an.tpr
-rw-r--r-- 1 romashka staff 5,4K 27 мая 19:52 et_an.tpr
In [28]:
!ls -lh *.tpr
-rw-r--r-- 1 romashka staff 5,4K 27 мая 19:52 et_an.tpr -rw-r--r-- 1 romashka staff 5,4K 27 мая 19:47 et_be.tpr -rw-r--r-- 1 romashka staff 5,4K 27 мая 19:50 et_nh.tpr -rw-r--r-- 1 romashka staff 5,4K 27 мая 19:50 et_sd.tpr -rw-r--r-- 1 romashka staff 5,4K 27 мая 19:50 et_vr.tpr
In [29]:
methods = ["vr", "nh", "an", "sd"]
for m in methods:
print("=" * 60)
print("Запускаем:", m)
!gmx mdrun -deffnm et_{m} -v -nt 1
============================================================
Запускаем: vr
:-) GROMACS - gmx mdrun, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx mdrun -deffnm et_vr -v -nt 1
Reading file et_vr.tpr, VERSION 2026.0-conda_forge (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.
Can not increase nstlist because verlet-buffer-tolerance is not set or used
Using 1 MPI thread
Using 1 OpenMP thread
starting mdrun 'one ethane molecule'
50000 steps, 50.0 ps.
step 50000, remaining wall clock time: 0 s
Writing final coordinates.
Core t (s) Wall t (s) (%)
Time: 0.140 0.140 100.0
(ns/day) (hour/ns) (ms/step) (Matom*steps/s)
Performance: 30826.082 0.001 0.003 2.854
GROMACS reminds you: "It has been discovered that C++ provides a remarkable facility for concealing the trivial details of a program - such as where its bugs are." (David Keppel)
============================================================
Запускаем: nh
:-) GROMACS - gmx mdrun, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx mdrun -deffnm et_nh -v -nt 1
Reading file et_nh.tpr, VERSION 2026.0-conda_forge (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.
Can not increase nstlist because verlet-buffer-tolerance is not set or used
Using 1 MPI thread
Using 1 OpenMP thread
starting mdrun 'one ethane molecule'
50000 steps, 50.0 ps.
step 50000, remaining wall clock time: 0 s
Writing final coordinates.
Core t (s) Wall t (s) (%)
Time: 0.159 0.159 100.0
(ns/day) (hour/ns) (ms/step) (Matom*steps/s)
Performance: 27191.407 0.001 0.003 2.518
GROMACS reminds you: "It has been discovered that C++ provides a remarkable facility for concealing the trivial details of a program - such as where its bugs are." (David Keppel)
============================================================
Запускаем: an
:-) GROMACS - gmx mdrun, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx mdrun -deffnm et_an -v -nt 1
Reading file et_an.tpr, VERSION 2026.0-conda_forge (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.
Can not increase nstlist because verlet-buffer-tolerance is not set or used
Using 1 MPI thread
Using 1 OpenMP thread
starting mdrun 'one ethane molecule'
50000 steps, 50.0 ps.
step 50000, remaining wall clock time: 0 s
Writing final coordinates.
Core t (s) Wall t (s) (%)
Time: 0.145 0.145 100.0
(ns/day) (hour/ns) (ms/step) (Matom*steps/s)
Performance: 29726.057 0.001 0.003 2.752
GROMACS reminds you: "I invented the term 'Object-Oriented', and I can tell you I did not have C++ in mind." (Alan Kay, author of Smalltalk)
============================================================
Запускаем: sd
:-) GROMACS - gmx mdrun, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx mdrun -deffnm et_sd -v -nt 1
Reading file et_sd.tpr, VERSION 2026.0-conda_forge (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.
Can not increase nstlist because verlet-buffer-tolerance is not set or used
Using 1 MPI thread
Using 1 OpenMP thread
starting mdrun 'one ethane molecule'
50000 steps, 50.0 ps.
step 49900, remaining wall clock time: 0 s
Writing final coordinates.
step 50000, remaining wall clock time: 0 s
Core t (s) Wall t (s) (%)
Time: 0.119 0.119 100.0
(ns/day) (hour/ns) (ms/step) (Matom*steps/s)
Performance: 36402.083 0.001 0.002 3.371
GROMACS reminds you: "I invented the term 'Object-Oriented', and I can tell you I did not have C++ in mind." (Alan Kay, author of Smalltalk)
In [30]:
!ls -lh et_*.trr et_*.edr et_*.log et_*.gro
-rw-r--r-- 1 romashka staff 978K 27 мая 19:52 et_an.edr -rw-r--r-- 1 romashka staff 609B 27 мая 19:52 et_an.gro -rw-r--r-- 1 romashka staff 241K 27 мая 19:52 et_an.log -rw-r--r-- 1 romashka staff 1,5M 27 мая 19:52 et_an.trr -rw-r--r-- 1 romashka staff 1,0M 27 мая 19:48 et_be.edr -rw-r--r-- 1 romashka staff 609B 27 мая 19:48 et_be.gro -rw-r--r-- 1 romashka staff 255K 27 мая 19:48 et_be.log -rw-r--r-- 1 romashka staff 1,5M 27 мая 19:48 et_be.trr -rw-r--r-- 1 romashka staff 403B 27 мая 19:44 et_box.gro -rw-r--r-- 1 romashka staff 997K 27 мая 19:52 et_nh.edr -rw-r--r-- 1 romashka staff 609B 27 мая 19:52 et_nh.gro -rw-r--r-- 1 romashka staff 255K 27 мая 19:52 et_nh.log -rw-r--r-- 1 romashka staff 1,5M 27 мая 19:52 et_nh.trr -rw-r--r-- 1 romashka staff 978K 27 мая 19:52 et_sd.edr -rw-r--r-- 1 romashka staff 609B 27 мая 19:52 et_sd.gro -rw-r--r-- 1 romashka staff 242K 27 мая 19:52 et_sd.log -rw-r--r-- 1 romashka staff 1,5M 27 мая 19:52 et_sd.trr -rw-r--r-- 1 romashka staff 1,0M 27 мая 19:52 et_vr.edr -rw-r--r-- 1 romashka staff 609B 27 мая 19:52 et_vr.gro -rw-r--r-- 1 romashka staff 255K 27 мая 19:52 et_vr.log -rw-r--r-- 1 romashka staff 1,5M 27 мая 19:52 et_vr.trr
визуализация и графики энергий
In [34]:
methods = ["be", "vr", "nh", "an", "sd"]
for m in methods:
print("=" * 60)
print("Делаю PDB для PyMOL:", m)
!printf "0\n" | gmx trjconv -f et_{m}.trr -s et_{m}.tpr -o et_{m}_vis.pdb -dt 1
============================================================
Делаю PDB для PyMOL: be
:-) GROMACS - gmx trjconv, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx trjconv -f et_be.trr -s et_be.tpr -o et_be_vis.pdb -dt 1
Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write pdb: Protein data bank file
Reading file et_be.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_be.tpr, VERSION 2026.0-conda_forge (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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000 -> frame 40 time 40.000
-> frame 50 time 50.000
Last written: frame 50 time 50.000
GROMACS reminds you: "Men love to wonder, and that is the seed of science." (Ralph Waldo Emerson)
============================================================
Делаю PDB для PyMOL: vr
:-) GROMACS - gmx trjconv, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx trjconv -f et_vr.trr -s et_vr.tpr -o et_vr_vis.pdb -dt 1
Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write pdb: Protein data bank file
Reading file et_vr.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_vr.tpr, VERSION 2026.0-conda_forge (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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000 -> frame 40 time 40.000
-> frame 50 time 50.000
Last written: frame 50 time 50.000
GROMACS reminds you: "Men love to wonder, and that is the seed of science." (Ralph Waldo Emerson)
============================================================
Делаю PDB для PyMOL: nh
:-) GROMACS - gmx trjconv, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx trjconv -f et_nh.trr -s et_nh.tpr -o et_nh_vis.pdb -dt 1
Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write pdb: Protein data bank file
Reading file et_nh.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_nh.tpr, VERSION 2026.0-conda_forge (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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000 -> frame 40 time 40.000
-> frame 50 time 50.000
Last written: frame 50 time 50.000
GROMACS reminds you: "Men love to wonder, and that is the seed of science." (Ralph Waldo Emerson)
============================================================
Делаю PDB для PyMOL: an
:-) GROMACS - gmx trjconv, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx trjconv -f et_an.trr -s et_an.tpr -o et_an_vis.pdb -dt 1
Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write pdb: Protein data bank file
Reading file et_an.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_an.tpr, VERSION 2026.0-conda_forge (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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000 -> frame 40 time 40.000
-> frame 50 time 50.000
Last written: frame 50 time 50.000
GROMACS reminds you: "Men love to wonder, and that is the seed of science." (Ralph Waldo Emerson)
============================================================
Делаю PDB для PyMOL: sd
:-) GROMACS - gmx trjconv, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx trjconv -f et_sd.trr -s et_sd.tpr -o et_sd_vis.pdb -dt 1
Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write pdb: Protein data bank file
Reading file et_sd.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_sd.tpr, VERSION 2026.0-conda_forge (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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000 -> frame 40 time 40.000
-> frame 50 time 50.000
Last written: frame 50 time 50.000
GROMACS reminds you: "Men love to wonder, and that is the seed of science." (Ralph Waldo Emerson)
In [36]:
!ls -lh *_vis.pdb
-rw-r--r-- 1 romashka staff 42K 27 мая 20:00 et_an_vis.pdb -rw-r--r-- 1 romashka staff 42K 27 мая 20:00 et_be_vis.pdb -rw-r--r-- 1 romashka staff 42K 27 мая 20:00 et_nh_vis.pdb -rw-r--r-- 1 romashka staff 42K 27 мая 20:00 et_sd_vis.pdb -rw-r--r-- 1 romashka staff 42K 27 мая 20:00 et_vr_vis.pdb
In [37]:
import pymol2
print("PyMOL импортировался нормально")
PyMOL импортировался нормально
In [38]:
import os
import pymol2
from IPython.display import Image, display
methods = ["be", "vr", "nh", "an", "sd"]
names = {
"be": "Berendsen",
"vr": "Velocity rescale",
"nh": "Nose-Hoover",
"an": "Andersen",
"sd": "Stochastic dynamics"
}
png_files = []
for m in methods:
pdb_file = f"et_{m}_vis.pdb"
png_file = f"et_{m}_pymol.png"
with pymol2.PyMOL() as pymol:
cmd = pymol.cmd
cmd.load(pdb_file, "ethane")
cmd.hide("everything")
cmd.show("sticks", "ethane")
cmd.show("spheres", "ethane")
cmd.set("stick_radius", 0.12)
cmd.set("sphere_scale", 0.25)
cmd.set("ray_opaque_background", 0)
cmd.bg_color("white")
cmd.orient("ethane")
cmd.zoom("ethane", 2.0)
cmd.png(png_file, width=900, height=700, dpi=200, ray=1)
png_files.append(png_file)
print(f"Сохранено: {png_file}")
Сохранено: et_be_pymol.png Сохранено: et_vr_pymol.png Сохранено: et_nh_pymol.png Сохранено: et_an_pymol.png Сохранено: et_sd_pymol.png
In [39]:
for m, png_file in zip(methods, png_files):
print(names[m])
display(Image(filename=png_file))
Berendsen
Velocity rescale
Nose-Hoover
Andersen
Stochastic dynamics
Для всех пяти методов контроля температуры траектории были конвертированы в формат PDB и визуализированы в PyMOL. Во всех случаях молекула этана сохраняет нормальную структуру: два атома углерода остаются связанными между собой, атомы водорода не отрываются, молекула не разрушается и не «разлетается». Это показывает, что топология была задана корректно, а выбранные параметры моделирования позволяют устойчиво считать систему.
На изображениях видно, что этан в разных методах находится в разных ориентациях. Это не означает принципиального изменения структуры молекулы, а отражает обычное поступательное/вращательное движение и тепловые флуктуации в ходе молекулярной динамики. По статичным картинкам нельзя строго судить о качестве термостата, но можно сделать предварительный вывод, что все методы поддерживают систему в физически разумном состоянии.
Визуально наиболее заметных нарушений структуры не наблюдается ни для Berendsen, ни для Velocity rescale, ни для Nose–Hoover, Andersen и stochastic dynamics. Основные различия между методами нужно дальше оценивать количественно — по графикам кинетической и потенциальной энергии, а также по распределению длины связи C–C. Именно эти параметры лучше покажут, какой термостат реалистичнее поддерживает температуру.
Анализ потенциальной и кинетической энергии
In [40]:
methods = ["be", "vr", "nh", "an", "sd"]
names = {
"be": "Berendsen",
"vr": "Velocity rescale",
"nh": "Nose-Hoover",
"an": "Andersen",
"sd": "Stochastic dynamics"
}
for m in methods:
print("=" * 60)
print("Извлекаю энергии для:", names[m])
!printf "Potential\nKinetic-En.\n0\n" | gmx energy -f et_{m}.edr -o et_{m}_en.xvg -xvg none
============================================================
Извлекаю энергии для: Berendsen
:-) GROMACS - gmx energy, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx energy -f et_be.edr -o et_be_en.xvg -xvg none
Opened et_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
Last energy frame read 5000 time 50.000
Statistics over 50001 steps [ 0.0000 through 50.0000 ps ], 2 data sets
All statistics are over 5001 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 20.6367 1.2 4.69968 -8.07328 (kJ/mol)
Kinetic En. 30.9125 0.79 4.63576 -0.708132 (kJ/mol)
GROMACS reminds you: "The best model of a cat is another cat..., specially the same cat." (Arturo Rosenblueth)
============================================================
Извлекаю энергии для: Velocity rescale
:-) GROMACS - gmx energy, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx energy -f et_vr.edr -o et_vr_en.xvg -xvg none
Opened et_vr.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
Last energy frame read 5000 time 50.000
Statistics over 50001 steps [ 0.0000 through 50.0000 ps ], 2 data sets
All statistics are over 5001 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 28.9386 0.39 6.64883 2.40135 (kJ/mol)
Kinetic En. 27.6476 0.28 8.3466 -0.168492 (kJ/mol)
GROMACS reminds you: "The best model of a cat is another cat..., specially the same cat." (Arturo Rosenblueth)
============================================================
Извлекаю энергии для: Nose-Hoover
:-) GROMACS - gmx energy, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx energy -f et_nh.edr -o et_nh_en.xvg -xvg none
Opened et_nh.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
Last energy frame read 5000 time 50.000
Statistics over 50001 steps [ 0.0000 through 50.0000 ps ], 2 data sets
All statistics are over 5001 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 27.8489 0.098 6.41751 0.347796 (kJ/mol)
Kinetic En. 28.8989 0.47 8.77194 1.8962 (kJ/mol)
GROMACS reminds you: "The best model of a cat is another cat..., specially the same cat." (Arturo Rosenblueth)
============================================================
Извлекаю энергии для: Andersen
:-) GROMACS - gmx energy, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx energy -f et_an.edr -o et_an_en.xvg -xvg none
Opened et_an.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 Temperature 12 Pressure
13 Vir-XX 14 Vir-XY 15 Vir-XZ 16 Vir-YX
17 Vir-YY 18 Vir-YZ 19 Vir-ZX 20 Vir-ZY
21 Vir-ZZ 22 Pres-XX 23 Pres-XY 24 Pres-XZ
25 Pres-YX 26 Pres-YY 27 Pres-YZ 28 Pres-ZX
29 Pres-ZY 30 Pres-ZZ 31 #Surf*SurfTen 32 T-System
Last energy frame read 5000 time 50.000
Statistics over 50001 steps [ 0.0000 through 50.0000 ps ], 2 data sets
All statistics are over 5001 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 30.8877 0.58 7.47101 -1.63671 (kJ/mol)
Kinetic En. 25.2858 0.73 7.77883 -1.85944 (kJ/mol)
GROMACS reminds you: "The best model of a cat is another cat..., specially the same cat." (Arturo Rosenblueth)
============================================================
Извлекаю энергии для: Stochastic dynamics
:-) GROMACS - gmx energy, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx energy -f et_sd.edr -o et_sd_en.xvg -xvg none
Opened et_sd.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 Temperature 12 Pressure
13 Vir-XX 14 Vir-XY 15 Vir-XZ 16 Vir-YX
17 Vir-YY 18 Vir-YZ 19 Vir-ZX 20 Vir-ZY
21 Vir-ZZ 22 Pres-XX 23 Pres-XY 24 Pres-XZ
25 Pres-YX 26 Pres-YY 27 Pres-YZ 28 Pres-ZX
29 Pres-ZY 30 Pres-ZZ 31 #Surf*SurfTen 32 T-System
Last energy frame read 5000 time 50.000
Statistics over 50001 steps [ 0.0000 through 50.0000 ps ], 2 data sets
All statistics are over 5001 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 32.171 0.2 7.82533 0.485701 (kJ/mol)
Kinetic En. 30.1759 0.19 8.44048 0.8254 (kJ/mol)
GROMACS reminds you: "The best model of a cat is another cat..., specially the same cat." (Arturo Rosenblueth)
In [41]:
!ls -lh *_en.xvg
-rw-r--r-- 1 romashka staff 181K 27 мая 20:05 et_an_en.xvg -rw-r--r-- 1 romashka staff 181K 27 мая 20:05 et_be_en.xvg -rw-r--r-- 1 romashka staff 181K 27 мая 20:05 et_nh_en.xvg -rw-r--r-- 1 romashka staff 181K 27 мая 20:05 et_sd_en.xvg -rw-r--r-- 1 romashka staff 181K 27 мая 20:05 et_vr_en.xvg
In [42]:
!head et_be_en.xvg
0.000000 9.788048 28.761471
0.010000 24.819202 13.817494
0.020000 19.229340 19.697229
0.030000 32.306797 8.498548
0.040000 17.574211 22.856560
0.050000 21.245382 19.915182
0.060000 14.364448 26.751240
0.070000 28.517815 13.300981
0.080000 23.093742 18.927780
0.090000 30.446293 13.218240
In [43]:
import numpy as np
import matplotlib.pyplot as plt
methods = ["be", "vr", "nh", "an", "sd"]
names = {
"be": "Berendsen",
"vr": "Velocity rescale",
"nh": "Nose-Hoover",
"an": "Andersen",
"sd": "Stochastic dynamics"
}
for m in methods:
data = np.loadtxt(f"et_{m}_en.xvg")
t = data[:, 0]
potential = data[:, 1]
kinetic = data[:, 2]
plt.figure(figsize=(8, 4))
plt.plot(t, potential, label="Potential energy")
plt.plot(t, kinetic, label="Kinetic energy")
plt.xlabel("Time, ps")
plt.ylabel("Energy, kJ/mol")
plt.title(f"Potential and kinetic energy: {names[m]}")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(f"energy_{m}.png", dpi=300)
plt.show()
Для всех пяти термостатов были построены графики потенциальной и кинетической энергии. Во всех случаях система остаётся стабильной: энергия не уходит бесконечно вверх или вниз, молекула этана не разрушается.
Метод Berendsen даёт более сглаженное поведение энергии, что связано с искусственным масштабированием скоростей. Поэтому он хорошо быстро подгоняет температуру, но хуже воспроизводит реальные тепловые флуктуации.
Методы Velocity rescale и Nose–Hoover дают более выраженные флуктуации энергии вокруг среднего значения, что выглядит физически более реалистично. Andersen и stochastic dynamics показывают более резкие колебания, так как в них есть случайная компонента.
Вывод: все методы поддерживают температуру, но Berendsen выглядит менее реалистичным из-за подавления флуктуаций. Более корректными по характеру поведения энергии выглядят Velocity rescale и Nose–Hoover.
Анализ распределения длины связи C–C
In [44]:
with open("b.ndx", "w") as f:
f.write("[ b ]\n")
f.write("1 2\n")
!cat b.ndx
[ b ] 1 2
In [45]:
methods = ["be", "vr", "nh", "an", "sd"]
for m in methods:
print("=" * 60)
print("Считаю длину связи C-C для:", m)
!gmx distance -f et_{m}.trr -s et_{m}.tpr -n b.ndx -select 'atomnr 1 plus atomnr 2' -oall dist_{m}.xvg -oh bond_{m}.xvg -xvg none
============================================================
Считаю длину связи C-C для: be
:-) GROMACS - gmx distance, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx distance -f et_be.trr -s et_be.tpr -n b.ndx -select 'atomnr 1 plus atomnr 2' -oall dist_be.xvg -oh bond_be.xvg -xvg none
NOTE: You provided an index file
b.ndx
(with -n), but it was not used by any selection.
Reading file et_be.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_be.tpr, VERSION 2026.0-conda_forge (single precision)
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000
Analyzed 5001 frames, last time 50.000
atomnr 1 plus atomnr 2:
Number of samples: 5001
Average distance: 0.15372 nm
Standard deviation: 0.00164 nm
GROMACS reminds you: "Der Ball ist rund, das Spiel dauert 90 minuten, alles andere ist Theorie" (Lola rennt)
============================================================
Считаю длину связи C-C для: vr
:-) GROMACS - gmx distance, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx distance -f et_vr.trr -s et_vr.tpr -n b.ndx -select 'atomnr 1 plus atomnr 2' -oall dist_vr.xvg -oh bond_vr.xvg -xvg none
NOTE: You provided an index file
b.ndx
(with -n), but it was not used by any selection.
Reading file et_vr.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_vr.tpr, VERSION 2026.0-conda_forge (single precision)
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000
Analyzed 5001 frames, last time 50.000
atomnr 1 plus atomnr 2:
Number of samples: 5001
Average distance: 0.15344 nm
Standard deviation: 0.00262 nm
GROMACS reminds you: "Der Ball ist rund, das Spiel dauert 90 minuten, alles andere ist Theorie" (Lola rennt)
============================================================
Считаю длину связи C-C для: nh
:-) GROMACS - gmx distance, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx distance -f et_nh.trr -s et_nh.tpr -n b.ndx -select 'atomnr 1 plus atomnr 2' -oall dist_nh.xvg -oh bond_nh.xvg -xvg none
NOTE: You provided an index file
b.ndx
(with -n), but it was not used by any selection.
Reading file et_nh.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_nh.tpr, VERSION 2026.0-conda_forge (single precision)
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000
Analyzed 5001 frames, last time 50.000
atomnr 1 plus atomnr 2:
Number of samples: 5001
Average distance: 0.15348 nm
Standard deviation: 0.00175 nm
GROMACS reminds you: "Der Ball ist rund, das Spiel dauert 90 minuten, alles andere ist Theorie" (Lola rennt)
============================================================
Считаю длину связи C-C для: an
:-) GROMACS - gmx distance, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx distance -f et_an.trr -s et_an.tpr -n b.ndx -select 'atomnr 1 plus atomnr 2' -oall dist_an.xvg -oh bond_an.xvg -xvg none
NOTE: You provided an index file
b.ndx
(with -n), but it was not used by any selection.
Reading file et_an.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_an.tpr, VERSION 2026.0-conda_forge (single precision)
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000
Analyzed 5001 frames, last time 50.000
atomnr 1 plus atomnr 2:
Number of samples: 5001
Average distance: 0.15331 nm
Standard deviation: 0.00324 nm
GROMACS reminds you: "Der Ball ist rund, das Spiel dauert 90 minuten, alles andere ist Theorie" (Lola rennt)
============================================================
Считаю длину связи C-C для: sd
:-) GROMACS - gmx distance, 2026.0-conda_forge (-:
Executable: /Users/romashka/Documents/jupyter/envs/jlab/bin.ARM_NEON_ASIMD/gmx
Data prefix: /Users/romashka/Documents/jupyter/envs/jlab
Working dir: /Users/romashka/Documents/jupyter/notebooks/task7_ethane_thermostats
Command line:
gmx distance -f et_sd.trr -s et_sd.tpr -n b.ndx -select 'atomnr 1 plus atomnr 2' -oall dist_sd.xvg -oh bond_sd.xvg -xvg none
NOTE: You provided an index file
b.ndx
(with -n), but it was not used by any selection.
Reading file et_sd.tpr, VERSION 2026.0-conda_forge (single precision)
Reading file et_sd.tpr, VERSION 2026.0-conda_forge (single precision)
trr version: GMX_trn_file (single precision)
Last frame 5000 time 50.000
Analyzed 5001 frames, last time 50.000
atomnr 1 plus atomnr 2:
Number of samples: 5001
Average distance: 0.15337 nm
Standard deviation: 0.00334 nm
GROMACS reminds you: "Der Ball ist rund, das Spiel dauert 90 minuten, alles andere ist Theorie" (Lola rennt)
In [46]:
!ls -lh dist_*.xvg bond_*.xvg
-rw-r--r-- 1 romashka staff 4,1K 27 мая 20:08 bond_an.xvg -rw-r--r-- 1 romashka staff 4,1K 27 мая 20:08 bond_be.xvg -rw-r--r-- 1 romashka staff 4,1K 27 мая 20:08 bond_nh.xvg -rw-r--r-- 1 romashka staff 4,1K 27 мая 20:08 bond_sd.xvg -rw-r--r-- 1 romashka staff 4,1K 27 мая 20:08 bond_vr.xvg -rw-r--r-- 1 romashka staff 103K 27 мая 20:08 dist_an.xvg -rw-r--r-- 1 romashka staff 103K 27 мая 20:08 dist_be.xvg -rw-r--r-- 1 romashka staff 103K 27 мая 20:08 dist_nh.xvg -rw-r--r-- 1 romashka staff 103K 27 мая 20:08 dist_sd.xvg -rw-r--r-- 1 romashka staff 103K 27 мая 20:08 dist_vr.xvg
In [47]:
import numpy as np
import matplotlib.pyplot as plt
methods = ["be", "vr", "nh", "an", "sd"]
names = {
"be": "Berendsen",
"vr": "Velocity rescale",
"nh": "Nose-Hoover",
"an": "Andersen",
"sd": "Stochastic dynamics"
}
for m in methods:
data = np.loadtxt(f"dist_{m}.xvg")
t = data[:, 0]
d = data[:, 1]
plt.figure(figsize=(6, 4))
plt.hist(d, bins=40)
plt.xlabel("C-C distance, nm")
plt.ylabel("Count")
plt.title(f"C-C bond length distribution: {names[m]}")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(f"bond_hist_{m}.png", dpi=300)
plt.show()
In [48]:
for m in methods:
data = np.loadtxt(f"dist_{m}.xvg")
d = data[:, 1]
print(names[m])
print("mean =", round(d.mean(), 4), "nm")
print("std =", round(d.std(), 4), "nm")
print()
Berendsen mean = 0.1537 nm std = 0.0017 nm Velocity rescale mean = 0.1534 nm std = 0.0026 nm Nose-Hoover mean = 0.1535 nm std = 0.0018 nm Andersen mean = 0.1533 nm std = 0.0033 nm Stochastic dynamics mean = 0.1534 nm std = 0.0034 nm
Для всех методов средняя длина связи C–C получилась почти одинаковой — около 0.153–0.154 нм, то есть молекула этана сохраняет нормальную структуру, а связь между атомами углерода не разрушается.
Различия видны по ширине распределения. Самое узкое распределение получилось для Berendsen и Nose–Hoover: стандартное отклонение около 0.0017–0.0018 нм. У Velocity rescale разброс больше — 0.0026 нм. Самые широкие распределения получились для Andersen и stochastic dynamics — около 0.0033–0.0034 нм, что связано со случайной компонентой этих методов.
Вывод: все термостаты сохраняют связь C–C, но Berendsen сильнее ограничивает флуктуации, а Andersen и stochastic dynamics дают более широкие колебания длины связи. Это согласуется с тем, что стохастические методы сильнее вмешиваются в движение атомов.
Сравнение быстродействия методов
In [50]:
import re
import pandas as pd
methods = ["be", "vr", "nh", "an", "sd"]
names = {
"be": "Berendsen",
"vr": "Velocity rescale",
"nh": "Nose-Hoover",
"an": "Andersen",
"sd": "Stochastic dynamics"
}
rows = []
for m in methods:
log_file = f"et_{m}.log"
with open(log_file) as f:
text = f.read()
# ищем строку Performance
match = re.search(r"Performance:\s+([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)", text)
if match:
ns_day = float(match.group(1))
hour_ns = float(match.group(2))
ms_step = float(match.group(3))
matom_steps_s = float(match.group(4))
else:
ns_day = None
hour_ns = None
ms_step = None
matom_steps_s = None
rows.append({
"method": names[m],
"ns/day": ns_day,
"hour/ns": hour_ns,
"ms/step": ms_step,
"Matom*steps/s": matom_steps_s
})
speed_table = pd.DataFrame(rows)
speed_table
Out[50]:
| method | ns/day | hour/ns | ms/step | Matom*steps/s | |
|---|---|---|---|---|---|
| 0 | Berendsen | 31286.161 | 0.001 | 0.003 | 2.897 |
| 1 | Velocity rescale | 30826.082 | 0.001 | 0.003 | 2.854 |
| 2 | Nose-Hoover | 27191.407 | 0.001 | 0.003 | 2.518 |
| 3 | Andersen | 29726.057 | 0.001 | 0.003 | 2.752 |
| 4 | Stochastic dynamics | 36402.083 | 0.001 | 0.002 | 3.371 |
In [51]:
speed_table.to_csv("speed_table.csv", index=False)
speed_table
Out[51]:
| method | ns/day | hour/ns | ms/step | Matom*steps/s | |
|---|---|---|---|---|---|
| 0 | Berendsen | 31286.161 | 0.001 | 0.003 | 2.897 |
| 1 | Velocity rescale | 30826.082 | 0.001 | 0.003 | 2.854 |
| 2 | Nose-Hoover | 27191.407 | 0.001 | 0.003 | 2.518 |
| 3 | Andersen | 29726.057 | 0.001 | 0.003 | 2.752 |
| 4 | Stochastic dynamics | 36402.083 | 0.001 | 0.002 | 3.371 |
In [52]:
import matplotlib.pyplot as plt
plt.figure(figsize=(8, 4))
plt.bar(speed_table["method"], speed_table["ns/day"])
plt.ylabel("Performance, ns/day")
plt.title("Simulation speed for different thermostats")
plt.xticks(rotation=30, ha="right")
plt.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.savefig("speed_comparison.png", dpi=300)
plt.show()
Для каждого метода контроля температуры была оценена скорость расчёта по значению Performance из log-файлов GROMACS. Быстродействие выражалось в ns/day: чем больше это значение, тем быстрее выполнялась симуляция.
По полученным данным все методы показали достаточно высокую скорость, так как система очень маленькая и состоит только из одной молекулы этана. Самым быстрым оказался метод stochastic dynamics, а самым медленным — Nose–Hoover. Методы Berendsen, Velocity rescale и Andersen дали близкие значения быстродействия.
Вывод: различия в скорости есть, но для такой маленькой системы они не являются главным критерием выбора термостата. Важнее учитывать физическую корректность метода: например, Berendsen работает быстро, но хуже воспроизводит правильные флуктуации, тогда как Velocity rescale и Nose–Hoover более реалистично поддерживают температуру.
Итоговое сравнение термостатов и общий вывод
In [53]:
import numpy as np
import pandas as pd
import re
methods = ["be", "vr", "nh", "an", "sd"]
names = {
"be": "Berendsen",
"vr": "Velocity rescale",
"nh": "Nose-Hoover",
"an": "Andersen",
"sd": "Stochastic dynamics"
}
rows = []
for m in methods:
# длина связи C-C
d = np.loadtxt(f"dist_{m}.xvg")[:, 1]
# скорость из log-файла
with open(f"et_{m}.log") as f:
text = f.read()
match = re.search(r"Performance:\s+([\d\.]+)", text)
speed = float(match.group(1)) if match else None
rows.append({
"Method": names[m],
"Mean C-C distance, nm": round(d.mean(), 4),
"Std C-C distance, nm": round(d.std(), 4),
"Performance, ns/day": round(speed, 1)
})
summary = pd.DataFrame(rows)
summary
Out[53]:
| Method | Mean C-C distance, nm | Std C-C distance, nm | Performance, ns/day | |
|---|---|---|---|---|
| 0 | Berendsen | 0.1537 | 0.0017 | 31286.2 |
| 1 | Velocity rescale | 0.1534 | 0.0026 | 30826.1 |
| 2 | Nose-Hoover | 0.1535 | 0.0018 | 27191.4 |
| 3 | Andersen | 0.1533 | 0.0033 | 29726.1 |
| 4 | Stochastic dynamics | 0.1534 | 0.0034 | 36402.1 |
In [54]:
summary.to_csv("summary_thermostats.csv", index=False)
В работе были сравнены пять методов контроля температуры в молекулярной динамике на примере одной молекулы этана: Berendsen, Velocity rescale, Nose–Hoover, Andersen и stochastic dynamics. Во всех случаях симуляции прошли стабильно: молекула этана не разрушалась, связь C–C сохранялась, а средняя длина этой связи оставалась около 0.153–0.154 нм.
По графикам энергии видно, что все методы поддерживают систему в устойчивом состоянии, однако характер флуктуаций различается. Berendsen даёт более сглаженное поведение энергии и более узкое распределение длины связи, то есть сильнее подавляет естественные флуктуации. Поэтому он удобен для быстрой предварительной стабилизации системы, но менее реалистичен для анализа равновесной динамики.
Velocity rescale и Nose–Hoover выглядят более физически корректными, потому что энергия флуктуирует вокруг среднего значения без явного ухода системы в нестабильное состояние. Andersen и stochastic dynamics дают более широкие распределения длины связи и более резкие колебания энергии, что связано со случайной компонентой этих методов.
По быстродействию самым быстрым в данном запуске оказался stochastic dynamics, а самым медленным — Nose–Hoover. Но так как система очень маленькая, различия в скорости не стоит считать главным критерием выбора метода. Более важным является то, насколько реалистично термостат воспроизводит тепловые флуктуации.
Общий вывод: для данной системы наиболее разумными методами контроля температуры выглядят Velocity rescale и Nose–Hoover. Berendsen менее реалистичен из-за подавления флуктуаций, а Andersen и stochastic dynamics сильнее вмешиваются в движение атомов за счёт случайной компоненты.
In [ ]: