Skip to content

Commit 89301b4

Browse files
simongravellesimongravelle
committed
updated bulk
1 parent f731f6c commit 89301b4

1 file changed

Lines changed: 131 additions & 58 deletions

File tree

‎docs/sphinx/source/tutorials/tutorial1/bulk-solution.rst‎

Lines changed: 131 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -902,28 +902,97 @@ at equilibrium temperature and pressure, **npt.gro**, using:
902902
gmx grompp -f inputs/production.mdp -c npt.gro -o production -pp production -po production
903903
gmx mdrun -v -deffnm production
904904
905-
Radial distribution function
906-
----------------------------
907-
908-
After the simulation is completed, let us compute the radial
909-
distribution functions between :math:`\text{Na}^+` and
910-
:math:`\text{H}_2\text{O}`, :math:`\text{SO}_4^{2-}` and
911-
:math:`\text{H}_2\text{O}`, as well as in between
912-
:math:`\text{H}_2\text{O}` molecules. This can be done using
905+
Measurement
906+
===========
907+
908+
After completing the simulation, we proceed to compute the radial distribution functions (rdf):
909+
910+
.. math::
911+
912+
g(r) = \frac{V}{N_{\text{ref}} \rho} \frac{dN(r)}{dr},
913+
914+
where :math:`V` is the volume of the simulation box, :math:`N_{\text{ref}}` is
915+
the number of reference atoms, :math:`\rho` is the average number density of
916+
particles in the system, and :math:`\frac{dN(r)}{dr}` is the number of particles
917+
in a spherical shell of thickness :math:`dr` around a reference particle at
918+
a distance :math:`r`.
919+
920+
First, let us measure the rdf between :math:`\text{Na}^+`
921+
ions and :math:`\text{H}_2\text{O}` molecules, as well as between :math:`\text{SO}_4^{2-}`
922+
ions and :math:`\text{H}_2\text{O}`. This can be done using
913923
the ``gmx rdf`` command as follows:
914924

915-
.. code-block:: bash
925+
.. code-block:: bw
916926
917-
gmx rdf -f production.xtc -s production.tpr -o na-sol-rdf.xvg
927+
gmx rdf -f production.xtc -s production.tpr -o production-rdf-na-h2o.xvg
918928
919-
Selecting the sodium ions, and then the water. Repeat the same operation for
920-
the sulfate and water, and for the water and water. For the water-water
921-
RDF, it is better to exclude the intra-molecular contribution using
922-
the *-excl* option, as follows:
929+
Then select the sodium ions as *reference* by typing 3, the water
930+
as *selection* by typing 4, and press ``Ctrl+D``. The same can be done
931+
for :math:`\text{SO}_4^{2-}` ions by typing:
923932

924-
.. code-block:: bash
933+
.. code-block:: bw
934+
935+
${gmx} rdf -f production.xtc -s production.tpr -o production-rdf-so4-h2o.xvg
936+
937+
and then by typing 2 and 4.
938+
939+
The results show...
940+
941+
The main issue with the calculated rdf, is that it includes all the atoms from
942+
thr :math:`\text{H}_2\text{O}` molecules (including the hydrogen atoms) and all
943+
the atoms from the :math:`\text{SO}_4^{2-}`, leading to more peaks and dephts
944+
and a more difficult analysis. Rdfs would be easiers to interpret, if only the
945+
water oxygen atoms (with name ``OW1``) and :math:`\text{SO}_4^{2-}` ions
946+
sulfur atoms (with name ``S1``) where included in the analysis. As these groups were not
947+
included in the original group, we have to create them ourself.
948+
949+
To create groups, we can use the ``gmx make_ndx`` command as follow:
950+
951+
.. code-block:: bw
952+
953+
gmx make_ndx -f production.tpr << EOF
954+
a OW1
955+
a S1
956+
q
957+
EOF
958+
959+
And then type ``a OW1`` and press enter, ``a S1`` and press enter, and then
960+
press ``q`` to save and quit. This will create a file name **index.ndx** that
961+
contain two more groups (named OW1 and S1) alongside the default ones:
962+
963+
.. code-block:: bw
964+
965+
(...)
966+
3223 3224 3225 3226 3227 3228 3229 3230 3231 3232 3233 3234 3235 3236 3237
967+
3238 3239 3240 3241 3242
968+
[ non-Water ]
969+
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
970+
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
971+
31 32 33 34 35 36 37 38 39 40 41 42
972+
[ OW1 ]
973+
43 47 51 55 59 63 67 71 75 79 83 87 91 95 99
974+
103 107 111 115 119 123 127 131 135 139 143 147 151 155 159
975+
(...)
976+
3163 3167 3171 3175 3179 3183 3187 3191 3195 3199 3203 3207 3211 3215 3219
977+
3223 3227 3231 3235 3239
978+
[ S1 ]
979+
5 10 15 20 25 30
980+
981+
Then, let us rerun the ``gmx rdf`` command using the **index.ndx** file, and
982+
selecting the newly created groups:
983+
984+
.. code-block:: bw
985+
986+
gmx rdf -f production.xtc -s production.tpr -o production-rdf-na-OW1.xvg -n index.ndx
987+
988+
and select 3 and 7.
989+
990+
.. code-block:: bw
991+
992+
gmx rdf -f production.xtc -s production.tpr -o production-rdf-so4-OW1.xvg -n index.ndx
993+
994+
and select 8 and 7.
925995

926-
gmx rdf -f production.xtc -s production.tpr -o sol-sol-rdf.xvg -excl
927996

928997
.. figure:: ../figures/level1/bulk-solution/rdf-production-light.png
929998
:alt: Gromacs tutorial RDF radial distribution function
@@ -943,56 +1012,60 @@ the different species of the fluid. For instance, it can be seen that
9431012
there is a strong hydration layer of water around sodium ions at a typical
9441013
distance of :math:`2.4 ~ \text{Å}` from the center of the sodium ion.
9451014

946-
Mean square displacement
947-
========================
1015+
.. include:: ../../non-tutorials/accessfile.rst
9481016

949-
To probe the system dynamics, let us compute the mean square
950-
displacement for all 3 species. For the sulfate ion, type:
9511017

952-
.. code-block:: bash
9531018

954-
gmx msd -f production.xtc -s production.tpr -o so4-msd.xvg
9551019

956-
Select the :math:`\text{SO}_4^{2-}` ions (in my case, it is done by typing
957-
*2*), and then press *ctrl D*. Repeat the same operation for the sodium
958-
ions and for the water molecules.
1020+
..
1021+
Mean square displacement
1022+
========================
9591023
960-
The slope of the MSD in the limit of long times gives an estimate of the diffusion
961-
coefficient, following :math:`D = \text{MSD} / 2 d t`,
962-
where :math:`d = 3` is the dimension of the system. Here,
963-
I find a value of :math:`1.4 \mathrm{e}-5 ~ \text{cm}^2/\text{s}` for the
964-
diffusion coefficient of the sulfur ions.
965-
966-
Repeat the same for :math:`\text{Na}^+` and water.
967-
968-
For sodium, I find a value of :math:`1.6 \mathrm{e}-5 ~ \text{cm}^2/\text{s}`
969-
for the diffusion coefficient,
970-
and for water :math:`5.3 \mathrm{e}-5 ~ \text{cm}^2/\text{s}`.
971-
For comparison, the experimental diffusion coefficient of *pure* water at
972-
temperature :math:`T = 360~\text{K}`
973-
is :math:`7.3 \mathrm{e}-5 ~ \text{cm}^2/\text{s}` :cite:`simpson1958diffusion`.
974-
In the presence of ions, the diffusion coefficient of water is expected to
975-
be reduced.
976-
977-
.. figure:: ../figures/level1/bulk-solution/msd-production-light.png
978-
:alt: Gromacs tutorial : diffusion coefficient
979-
:class: only-light
1024+
To probe the system dynamics, let us compute the mean square
1025+
displacement for all 3 species. For the sulfate ion, type:
9801026

981-
.. figure:: ../figures/level1/bulk-solution/msd-production-dark.png
982-
:alt: Gromacs tutorial : diffusion coefficient
983-
:class: only-dark
1027+
.. code-block:: bash
9841028
985-
.. container:: figurelegend
1029+
gmx msd -f production.xtc -s production.tpr -o so4-msd.xvg
9861030
987-
Figure: Mean square displacement (msd) for the three species. The dashed line
988-
highlight the proportionality between msd and time :math:`t` which is expected
989-
at long times, when the system reaches the diffusive regime.
1031+
Select the :math:`\text{SO}_4^{2-}` ions (in my case, it is done by typing
1032+
*2*), and then press *ctrl D*. Repeat the same operation for the sodium
1033+
ions and for the water molecules.
9901034

991-
.. admonition:: About diffusion coefficient measurement in molecular simulations
992-
:class: info
1035+
The slope of the MSD in the limit of long times gives an estimate of the diffusion
1036+
coefficient, following :math:`D = \text{MSD} / 2 d t`,
1037+
where :math:`d = 3` is the dimension of the system. Here,
1038+
I find a value of :math:`1.4 \mathrm{e}-5 ~ \text{cm}^2/\text{s}` for the
1039+
diffusion coefficient of the sulfur ions.
1040+
1041+
Repeat the same for :math:`\text{Na}^+` and water.
9931042

994-
In principle, diffusion coefficients obtained from the MSD in a
995-
finite-sized box must be corrected, but this is beyond the scope of
996-
the present tutorial :cite:`loche2021transferable`.
1043+
For sodium, I find a value of :math:`1.6 \mathrm{e}-5 ~ \text{cm}^2/\text{s}`
1044+
for the diffusion coefficient,
1045+
and for water :math:`5.3 \mathrm{e}-5 ~ \text{cm}^2/\text{s}`.
1046+
For comparison, the experimental diffusion coefficient of *pure* water at
1047+
temperature :math:`T = 360~\text{K}`
1048+
is :math:`7.3 \mathrm{e}-5 ~ \text{cm}^2/\text{s}` :cite:`simpson1958diffusion`.
1049+
In the presence of ions, the diffusion coefficient of water is expected to
1050+
be reduced.
9971051

998-
.. include:: ../../non-tutorials/accessfile.rst
1052+
.. figure:: ../figures/level1/bulk-solution/msd-production-light.png
1053+
:alt: Gromacs tutorial : diffusion coefficient
1054+
:class: only-light
1055+
1056+
.. figure:: ../figures/level1/bulk-solution/msd-production-dark.png
1057+
:alt: Gromacs tutorial : diffusion coefficient
1058+
:class: only-dark
1059+
1060+
.. container:: figurelegend
1061+
1062+
Figure: Mean square displacement (msd) for the three species. The dashed line
1063+
highlight the proportionality between msd and time :math:`t` which is expected
1064+
at long times, when the system reaches the diffusive regime.
1065+
1066+
.. admonition:: About diffusion coefficient measurement in molecular simulations
1067+
:class: info
1068+
1069+
In principle, diffusion coefficients obtained from the MSD in a
1070+
finite-sized box must be corrected, but this is beyond the scope of
1071+
the present tutorial :cite:`loche2021transferable`.

0 commit comments

Comments
 (0)