In [1]:
# Сделаем файл etane.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 [2]:
# Построим файл топологии для этана et.top
In [ ]:
#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 6 1
# 1 2 7 1
# 1 2 8 1
#[ dihedrals ]
#; ai aj ak al funct
# 3 1 2 6 3
# 3 1 2 7 3
# 3 1 2 8 3
# 4 1 2 6 3
# 4 1 2 7 3
# 4 1 2 8 3
# 5 1 2 6 3
# 5 1 2 7 3
# 5 1 2 8 3
#[ pairs ]
#; список атомов 1-4
#; ai aj funct
# 3 6
# 3 7
# 3 8
# 4 6
# 4 7
# 4 8
# 5 6
# 5 7
# 5 8
#[ System ]
#; any text here
#first one
#[ molecules ]
#;Name count
# et 1
In [ ]:
#Вам даны 5 конфигурационных файлов с разными алгоритмами контроля температуры:
#be.mdp - метод Берендсена для контроля температуры.
#vr.mdp - метод "Velocity rescale" для контроля температуры.
#nh.mdp - метод Нуза-Хувера для контроля температуры.
#an.mdp - метод Андерсена для контроля температуры.
#sd.mdp - метод стохастической молекулярной динамики.
In [10]:
l = ['be', 'vr', 'nh', 'sd']
for s in l:
!gmx grompp -maxwarn 1 -f {s}.mdp -c etane.gro -p et.top -o et_{s}.tpr
:-) GROMACS - gmx grompp, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx grompp -maxwarn 1 -f be.mdp -c etane.gro -p et.top -o et_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, consider
setting nstcomm equal to nstcalcenergy for less overhead
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 -139038761
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 2 [file et.top, line 72]:
System has non-zero total charge: -0.309500
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 3 [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
NOTE 4 [file be.mdp]:
COM removal frequency is set to (1).
Other settings require a global communication frequency of 10.
Note that this will require additional global communication steps,
which will reduce performance when using multiple ranks.
Consider setting nstcomm to a multiple of 10.
There were 4 notes
There was 1 warning
Back Off! I just backed up et_be.tpr to ./#et_be.tpr.1#
GROMACS reminds you: "With four parameters I can fit an elephant, and with five I can make him wiggle his trunk." (John von Neumann)
:-) GROMACS - gmx grompp, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx grompp -maxwarn 1 -f vr.mdp -c etane.gro -p et.top -o et_vr.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 vr.mdp]:
nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, consider
setting nstcomm equal to nstcalcenergy for less overhead
Setting the LD random seed to -137512329
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 2 [file et.top, line 72]:
System has non-zero total charge: -0.309500
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 3 [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 10 Mb of data
NOTE 4 [file vr.mdp]:
COM removal frequency is set to (1).
Other settings require a global communication frequency of 10.
Note that this will require additional global communication steps,
which will reduce performance when using multiple ranks.
Consider setting nstcomm to a multiple of 10.
There were 4 notes
Back Off! I just backed up et_vr.tpr to ./#et_vr.tpr.2#
GROMACS reminds you: "When using an abacus, a human can achieve about 0.1 flops/watt. Super-computers achieve about 2 gigaflops/watt." (John Linford)
:-) GROMACS - gmx grompp, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx grompp -maxwarn 1 -f nh.mdp -c etane.gro -p et.top -o et_nh.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 nh.mdp]:
nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, consider
setting nstcomm equal to nstcalcenergy for less overhead
NOTE 2 [file nh.mdp]:
leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1
Setting the LD random seed to -319554593
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 et.top, line 72]:
System has non-zero total charge: -0.309500
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 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 10 Mb of data
NOTE 5 [file nh.mdp]:
COM removal frequency is set to (1).
Other settings require a global communication frequency of 10.
Note that this will require additional global communication steps,
which will reduce performance when using multiple ranks.
Consider setting nstcomm to a multiple of 10.
There were 5 notes
Back Off! I just backed up et_nh.tpr to ./#et_nh.tpr.2#
GROMACS reminds you: "When using an abacus, a human can achieve about 0.1 flops/watt. Super-computers achieve about 2 gigaflops/watt." (John Linford)
:-) GROMACS - gmx grompp, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx grompp -maxwarn 1 -f sd.mdp -c etane.gro -p et.top -o et_sd.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 sd.mdp]:
nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, consider
setting nstcomm equal to nstcalcenergy for less overhead
Setting the LD random seed to -545792009
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 2 [file et.top, line 72]:
System has non-zero total charge: -0.309500
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 3 [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 10 Mb of data
NOTE 4 [file sd.mdp]:
COM removal frequency is set to (1).
Other settings require a global communication frequency of 10.
Note that this will require additional global communication steps,
which will reduce performance when using multiple ranks.
Consider setting nstcomm to a multiple of 10.
There were 4 notes
Back Off! I just backed up et_sd.tpr to ./#et_sd.tpr.2#
GROMACS reminds you: "Science Won't Change You" (The Talking Heads)
In [11]:
# У нас должно получиться 4 tpr файлов. (an убрали, так как в новом GROMACS он больше не поддердживается) Теперь для каждого из них запустим mdrun.
In [12]:
l = ['be', 'vr', 'nh', 'sd']
for s in l:
!gmx mdrun -deffnm et_{s} -v -nt 1
:-) GROMACS - gmx mdrun, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx mdrun -deffnm et_be -v -nt 1
Compiled SIMD: SSE4.1, but for this host/run AVX2_256 might be better (see
log).
Reading file et_be.tpr, VERSION 2022.5-Debian_2022.5_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.747
Using 1 MPI thread
Using 1 OpenMP thread
NOTE: OS CPU limit is lower than logical cpu count, thread pinning disabled.
starting mdrun 'first one'
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: 2.738 2.738 100.0
(ns/day) (hour/ns)
Performance: 7889.160 0.003
GROMACS reminds you: "This Puke Stinks Like Beer" (LIVE)
:-) GROMACS - gmx mdrun, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx mdrun -deffnm et_vr -v -nt 1
Compiled SIMD: SSE4.1, but for this host/run AVX2_256 might be better (see
log).
Reading file et_vr.tpr, VERSION 2022.5-Debian_2022.5_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.747
Using 1 MPI thread
Using 1 OpenMP thread
NOTE: OS CPU limit is lower than logical cpu count, thread pinning disabled.
starting mdrun 'first one'
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: 2.665 2.665 100.0
(ns/day) (hour/ns)
Performance: 8104.661 0.003
GROMACS reminds you: "I couldn't give a shit about ribosomes." (Björn Forsberg, presenting his thesis, including two papers on ribosomes)
:-) GROMACS - gmx mdrun, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx mdrun -deffnm et_nh -v -nt 1
Compiled SIMD: SSE4.1, but for this host/run AVX2_256 might be better (see
log).
Reading file et_nh.tpr, VERSION 2022.5-Debian_2022.5_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.747
Using 1 MPI thread
Using 1 OpenMP thread
NOTE: OS CPU limit is lower than logical cpu count, thread pinning disabled.
starting mdrun 'first one'
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: 2.682 2.682 100.0
(ns/day) (hour/ns)
Performance: 8052.763 0.003
GROMACS reminds you: "Bailed Out Of Edge Synchronization After 10,000 Iterations" (X/Motif)
:-) GROMACS - gmx mdrun, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx mdrun -deffnm et_sd -v -nt 1
Compiled SIMD: SSE4.1, but for this host/run AVX2_256 might be better (see
log).
Reading file et_sd.tpr, VERSION 2022.5-Debian_2022.5_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.747
Using 1 MPI thread
Using 1 OpenMP thread
NOTE: OS CPU limit is lower than logical cpu count, thread pinning disabled.
starting mdrun 'first one'
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: 2.744 2.744 100.0
(ns/day) (hour/ns)
Performance: 7871.285 0.003
GROMACS reminds you: "Don't waste pure thoughts on dirty enzymes." (Efraim Racker)
In [13]:
#Теперь переходим к анализу результатов. Начнем с визуального анализа.
#Для каждой из 4 систем проведите конвертацию в pdb и просмотрите в PyMol.
In [15]:
l = ['be', 'vr', 'nh', 'sd']
for s in l:
!echo 0 | gmx trjconv -f et_{s}.trr -s et_{s}.tpr -o et_{s}.pdb
:-) GROMACS - gmx trjconv, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx trjconv -f et_be.trr -s et_be.tpr -o et_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 et_be.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
Reading file et_be.tpr, VERSION 2022.5-Debian_2022.5_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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Reading frame 200 time 200.000 -> frame 199 time 199.000
Last written: frame 250 time 250.000
GROMACS reminds you: "It's just the way this stuff is done" (Built to Spill)
:-) GROMACS - gmx trjconv, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx trjconv -f et_vr.trr -s et_vr.tpr -o et_vr.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 et_vr.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
Reading file et_vr.tpr, VERSION 2022.5-Debian_2022.5_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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Reading frame 200 time 200.000 -> frame 199 time 199.000
Last written: frame 250 time 250.000
GROMACS reminds you: "It's just the way this stuff is done" (Built to Spill)
:-) GROMACS - gmx trjconv, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx trjconv -f et_nh.trr -s et_nh.tpr -o et_nh.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 et_nh.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
Reading file et_nh.tpr, VERSION 2022.5-Debian_2022.5_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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Reading frame 200 time 200.000 -> frame 199 time 199.000
Last written: frame 250 time 250.000
GROMACS reminds you: "It's just the way this stuff is done" (Built to Spill)
:-) GROMACS - gmx trjconv, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx trjconv -f et_sd.trr -s et_sd.tpr -o et_sd.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 et_sd.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
Reading file et_sd.tpr, VERSION 2022.5-Debian_2022.5_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: Selected 0: 'System'
trr version: GMX_trn_file (single precision)
Reading frame 200 time 200.000 -> frame 199 time 199.000
Last written: frame 250 time 250.000
GROMACS reminds you: "The aim of science is not to open the door to infinite wisdom, but to set a limit to infinite error." (Bertolt Brecht, Life of Galileo)
In [25]:
from IPython.display import Image
In [26]:
Image(filename='be.gif')
Out[26]:
<IPython.core.display.Image object>
In [27]:
Image(filename='vr.gif')
Out[27]:
<IPython.core.display.Image object>
In [28]:
Image(filename='nh.gif')
Out[28]:
<IPython.core.display.Image object>
In [29]:
Image(filename='sd.gif')
Out[29]:
<IPython.core.display.Image object>
In [16]:
#Сравним потенциальную энергию и кинетическую энергию для каждой из 4 систем.
In [19]:
l = ['be', 'vr', 'nh', 'sd']
for s in l:
!echo "8 9 0" | gmx energy -f et_{s}.edr -o et_{s}_en.xvg -xvg none
:-) GROMACS - gmx energy, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
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
Back Off! I just backed up et_be_en.xvg to ./#et_be_en.xvg.1#
Last energy frame read 2500 time 250.000
Statistics over 250001 steps [ 0.0000 through 250.0000 ps ], 2 data sets
All statistics are over 2501 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 16.8354 0.12 1.63759 0.678422 (kJ/mol)
Kinetic En. 26.06 0.05 1.59831 -0.22682 (kJ/mol)
GROMACS reminds you: "Time is the best appraiser of scientific work, and I am aware that an industrial discovery rarely produces all its fruit in the hands of its first inventor." (Louis Pasteur)
:-) GROMACS - gmx energy, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
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
Back Off! I just backed up et_vr_en.xvg to ./#et_vr_en.xvg.1#
Last energy frame read 2500 time 250.000
Statistics over 250001 steps [ 0.0000 through 250.0000 ps ], 2 data sets
All statistics are over 2501 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 16.9514 0.4 8.70699 1.34038 (kJ/mol)
Kinetic En. 26.8626 0.24 8.47283 0.971287 (kJ/mol)
GROMACS reminds you: "Time is the best appraiser of scientific work, and I am aware that an industrial discovery rarely produces all its fruit in the hands of its first inventor." (Louis Pasteur)
:-) GROMACS - gmx energy, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
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
Back Off! I just backed up et_nh_en.xvg to ./#et_nh_en.xvg.1#
Last energy frame read 2500 time 250.000
Statistics over 250001 steps [ 0.0000 through 250.0000 ps ], 2 data sets
All statistics are over 2501 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 25.8151 1.1 83.0321 -4.24394 (kJ/mol)
Kinetic En. 32.3337 0.81 71.5938 -0.771759 (kJ/mol)
GROMACS reminds you: "Time is the best appraiser of scientific work, and I am aware that an industrial discovery rarely produces all its fruit in the hands of its first inventor." (Louis Pasteur)
:-) GROMACS - gmx energy, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
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
Back Off! I just backed up et_sd_en.xvg to ./#et_sd_en.xvg.1#
Last energy frame read 2500 time 250.000
Statistics over 250001 steps [ 0.0000 through 250.0000 ps ], 2 data sets
All statistics are over 2501 points (frames)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential 13.0362 0.37 7.43987 1.57089 (kJ/mol)
Kinetic En. 26.2428 0.22 7.80083 0.512805 (kJ/mol)
GROMACS reminds you: "Time is the best appraiser of scientific work, and I am aware that an industrial discovery rarely produces all its fruit in the hands of its first inventor." (Louis Pasteur)
In [20]:
#Постройте графики изменения энергий. Рекомендуемый вид это lines. Графики добавьте в отчёт. Удобно загружать данные с numpy.loadtxt
In [22]:
import numpy as np
import matplotlib.pyplot as plt
In [27]:
l = ['be', 'vr', 'nh', 'sd']
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()
for i, s in enumerate(l):
a = np.loadtxt(f'et_{s}_en.xvg')
t = a[:, 0]
pot = a[:, 1]
kin = a[:, 2]
axes[i].plot(t, pot, label='Потенциальная', color='blue')
axes[i].plot(t, kin, label='Кинетическая', color='red', alpha=0.7)
axes[i].set_title(f'Алгоритм: {s.upper()}')
axes[i].set_xlabel('Время (пикосекунды)')
axes[i].set_ylabel('Энергия (кДж/моль)')
axes[i].legend()
axes[i].grid(True)
plt.tight_layout()
plt.show()
In [1]:
#Рассмотрим распределение длинны связи С-С за время моделирования.
#Сначала создадим индекс файл с одной связью. Файл: b.ndx
#С содержимым:
#[ b ]
#1 2
In [2]:
#Запустим утилиту по анализу связей distance
In [6]:
l = ['be', 'vr', 'nh', 'sd']
for s in l:
!echo 0 | gmx distance -f et_{s}.trr -s et_{s}.tpr -oh bond_{s}.xvg -n b.ndx -xvg none
:-) GROMACS - gmx distance, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx distance -f et_be.trr -s et_be.tpr -oh bond_be.xvg -n b.ndx -xvg none
Reading file et_be.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
Reading file et_be.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
trr version: GMX_trn_file (single precision)
Reading frame 200 time 200.000
Analyzed 251 frames, last time 250.000
b:
Number of samples: 251
Average distance: 0.15316 nm
Standard deviation: 0.00170 nm
GROMACS reminds you: "An expert is a person who has made all the mistakes that can be made in a very narrow field." (Niels Bohr)
:-) GROMACS - gmx distance, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx distance -f et_vr.trr -s et_vr.tpr -oh bond_vr.xvg -n b.ndx -xvg none
Reading file et_vr.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
Reading file et_vr.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
trr version: GMX_trn_file (single precision)
Reading frame 200 time 200.000
Analyzed 251 frames, last time 250.000
b:
Number of samples: 251
Average distance: 0.15336 nm
Standard deviation: 0.00378 nm
GROMACS reminds you: "An expert is a person who has made all the mistakes that can be made in a very narrow field." (Niels Bohr)
:-) GROMACS - gmx distance, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx distance -f et_nh.trr -s et_nh.tpr -oh bond_nh.xvg -n b.ndx -xvg none
Reading file et_nh.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
Reading file et_nh.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
trr version: GMX_trn_file (single precision)
Reading frame 200 time 200.000
Analyzed 251 frames, last time 250.000
b:
Number of samples: 251
Average distance: 0.14812 nm
Standard deviation: 0.00697 nm
GROMACS reminds you: "An expert is a person who has made all the mistakes that can be made in a very narrow field." (Niels Bohr)
:-) GROMACS - gmx distance, 2022.5-Debian_2022.5_2 (-:
Executable: /bin/gmx
Data prefix: /usr
Working dir: /home/students/y22/urazov.d/term8/Pr6
Command line:
gmx distance -f et_sd.trr -s et_sd.tpr -oh bond_sd.xvg -n b.ndx -xvg none
Reading file et_sd.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
Reading file et_sd.tpr, VERSION 2022.5-Debian_2022.5_2 (single precision)
trr version: GMX_trn_file (single precision)
Reading frame 200 time 200.000
Analyzed 251 frames, last time 250.000
b:
Number of samples: 251
Average distance: 0.15263 nm
Standard deviation: 0.00331 nm
GROMACS reminds you: "An expert is a person who has made all the mistakes that can be made in a very narrow field." (Niels Bohr)
In [8]:
#Построим гистограммы распределения длинн связей
In [9]:
import pandas as pd
import matplotlib.pyplot as plt
In [14]:
def hist(filename):
data = pd.read_csv(filename, sep='\s+', comment='@', header=None, names=['Time', 'Distance'])
plt.figure(figsize=(8, 5))
plt.hist(data['Distance'], bins=50, color='skyblue', edgecolor='black', alpha=0.7)
plt.title(f'Распределение расстояний: {filename}')
plt.xlabel('Расстояние (нм)')
plt.ylabel('Частота')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.show()
In [16]:
hist('bond_be.xvg')
In [18]:
hist('bond_vr.xvg')
In [15]:
hist('bond_nh.xvg')
In [19]:
hist('bond_sd.xvg')
In [24]:
#Из полученные результатов видно, что алгоритмы различаются дисперсией расстояний между C-C.
#У be (термостат Берендсена) - очень большая дисперсия
# vr (Velocity Rescaling) - плохо соотносится с распределением Больцмана
# nh (термостат Нозе-Гувера) - хороший алгоритм, малая дисперсия, лучше соотносится с распределением Больцмана
# sd (стохастическая динамика) - из-за симулированных постоянных толчков sd может не подходить для некоторых задач, но идеально подходит для молекулярной динамики в растворителе
#Таблица зависимости быстродействия от алгоритма:
# Алгоритм Быстродействие(ns/day)
# be 7889.160
# vr 8104.661
# nh 8052.763
# sd 7871.285
In [ ]: