.. _bulk-solution-label:

Bulk salt solution
******************

.. container:: hatnote

    The very basics of GROMACS through a simple example

.. figure:: ../figures/level1/bulk-solution/video-solution-white.webp
    :alt: Water solution of SO\ :sub:`4`\ :sup:`2-` and Na\ :sup:`+` ions visualized with VMD
    :class: only-light
    :height: 250
    :align: right

.. figure:: ../figures/level1/bulk-solution/video-solution-dark.webp
    :alt: Water solution of SO\ :sub:`4`\ :sup:`2-` and Na\ :sup:`+` ions visualized with VMD
    :class: only-dark
    :height: 250
    :align: right

..  container:: justify

    The objective of this tutorial is to use the
    open-source code GROMACS to perform a simple molecular
    dynamics simulation. The system is a bulk solution of water mixed
    with sodium (Na\ :sup:`+`) and sulfate
    (SO\ :sub:`4`\ :sup:`2-`) ions. 
    
..  container:: justify

    This tutorial illustrates several major ingredients of molecular
    dynamics simulations, such as energy minimization,
    thermostating, NVT and NPT equilibration, and
    trajectory visualization.

.. include:: ../../non-tutorials/needhelp.rst
.. include:: ../../non-tutorials/GROMACS2024.2.rst

The input files
===============

..  container:: justify

    In order to run the present simulation using GROMACS,
    we need the 3 following files (or sets of files):

    - 1) A **configuration file** (.gro) containing the
         initial positions of the atoms and the box
         dimensions.
    - 2) A **topology file** (.top) specifying the
         location of the force field files (.itp).
    - 3) An **input file** (.mdp) containing the
         parameters of the simulation (e.g. temperature, timestep).

..  container:: justify

    The specificity of the present tutorial is that both configuration
    and topology files were prepared with homemade Python scripts, see
    :ref:`create-conf-label`. In principle, it is also possible to
    prepare the system using GROMACS functionalities, such as 
    *gmx pdb2gmx*, *gmx trjconv*, or *gmx solvate*. This will be done
    in the next tutorial, :ref:`protein_electrolyte-label`.

1) The configuration file (.gro)
--------------------------------

..  container:: justify

    For the present simulation, the initial atom
    positions and box size are given in a *conf.gro* file
    (Gromos87 format) that you can download by clicking |conf-SO4.gro|.
    Save the *conf.gro* file in a folder. 

.. |conf-SO4.gro| raw:: html

    <a href="https://raw.githubusercontent.com/gromacstutorials/gromacstutorials-inputs/main/level1/bulk-solution/conf.gro" target="_blank">here</a>

..  container:: justify

    The *conf.gro* file looks like this:

..  code-block:: bw

    Na2SO4 solution
    2846
       1  SO4   O1    1   2.608   3.089   2.389
       1  SO4   O2    2   2.562   3.181   2.150
       1  SO4   O3    3   2.388   3.217   2.339
       1  SO4   O4    4   2.425   2.980   2.241
       1  SO4   S1    5   2.496   3.117   2.280
    (...)
     719  Sol  OW1 2843   3.220   2.380   1.540
     719  Sol  HW1 2844   3.279   2.456   1.540
     719  Sol  HW2 2845   3.279   2.304   1.540
     719  Sol  MW1 2846   3.230   2.380   1.540
   3.36000   3.36000   3.36000

..  container:: justify

    The first line *Na2SO4 solution* is just a comment, the second line is the
    total number of atoms, and the last line is the box dimension in nanometer,
    here 3.36 nm by 3.36 nm by 3.36 nm. Between the second and the
    last lines, there is one line per
    atom. Each line indicates, from left to right:

    - the residue ID, with all the atoms from the same SO\ :sub:`4`\ :sup:`2-` ion sharing the same residue ID,
    - the residue name,
    - the atom name,
    - the atom ID,
    - the atom position: *x*, *y*, and *z* coordinates in nanometer.
    
..  container:: justify

    Note that the format of a *.gro* file is fixed,
    and each column is in a fixed position.

..  container:: justify

    The *conf.gro* file can be visualized using VMD by
    typing in the terminal:

..  code-block:: bash

     vmd conf.gro

.. figure:: ../figures/level1/bulk-solution/step0-light.png
    :alt: Gromacs initial configuration of SO\ :sub:`4`\ :sup:`2-` and Na\ :sup:`+` ions visualized with VMD
    :class: only-light

.. figure:: ../figures/level1/bulk-solution/step0-dark.png
    :alt: Gromacs initial configuration of SO\ :sub:`4`\ :sup:`2-` and Na\ :sup:`+` ions visualized with VMD
    :class: only-dark

.. container:: figurelegend

    Figure: SO\ :sub:`4`\ :sup:`2-` ions, Na\ :sup:`+` ions, and water molecules.
    Oxygen atoms are in red, hydrogen in white, sodium in blue, and sulfur in yellow. For better rendering, the atom representation and colors
    were modified with respect to the default VMD representation. If 
    you want to obtain the same rendering, you can follow this 
    |vmd-tutorial| from the LAMMPS tutorials webpage to obtain a similar rendering.

.. |vmd-tutorial| raw:: html

    <a href="https://lammpstutorials.github.io/sphinx/build/html/tutorials/vmd/vmd-tutorial.html" target="_blank">VMD tutorial</a>

..  container:: justify

    As can be seen using VMD, the water molecules are
    arranged in a quite unrealistic and regular manner, with
    all dipoles facing in the same direction, and possibly
    wrong distances between some of the molecules and ions.
    This will be fixed during energy minimization.

2) The topology files (.top .itp)
-------------------------------------

..  container:: justify

    The topology file contains information about the
    interactions of the different atoms and molecules. You
    can download it by clicking |topol-SO4.top|.
    Place it in the same folder as the *conf.gro* file. The
    *topol.top* file looks like that:

.. |topol-SO4.top| raw:: html

    <a href="https://raw.githubusercontent.com/gromacstutorials/gromacstutorials-inputs/main/level1/bulk-solution/topol.top" target="_blank">here</a>

..  code-block:: bw

    #include "ff/forcefield.itp"
    #include "ff/h2o.itp"
    #include "ff/na.itp"
    #include "ff/so4.itp"

    [ System ]
    Na2SO4 solution

    [ Molecules ]
    SO4 6
    Na 12
    SOL 701

..  container:: justify

    The 4 first lines are used to include the values of the
    parameters that are given in 4 separate files located in the *ff/*
    folder (see below). 
    
..  container:: justify

    The rest of the *topol.top* file contains the system
    name (*Na2SO4 solution*), and the list of the residues. Here there
    is 6 SO\ :sub:`4`\ :sup:`2-` ions,
    12 Na\ :sup:`+` ions,
    and 701 water molecules. It is crucial
    that the order and number of residues in the topology file 
    match the order of the *conf.gro* file. If you open the *conf.gro*
    file, you can see that indeed the first 30 lines beyond to the
    6 SO\ :sub:`4`\ :sup:`2-` residues, the next 12 lines to the
    12 Na\ :sup:`+` residues, and that the remaining lines concern the
    water molecules.

..  container:: justify
    
    Create a folder named *ff/* next to the *conf.gro* and the
    *topol.top* files, and copy
    |forcefield.itp|, |h2o.itp|, |na.itp|, and |so4.itp|
    in it. These four files contain information about
    the atoms (names, masses, changes, Lennard-Jones
    coefficients) and residues (bond and angular
    constraints) for all the species that are involved here.

.. |forcefield.itp| raw:: html

    <a href="https://raw.githubusercontent.com/gromacstutorials/gromacstutorials-inputs/main/level1/bulk-solution/ff/forcefield.itp" target="_blank">forcefield.itp</a>

.. |h2o.itp| raw:: html

    <a href="https://raw.githubusercontent.com/gromacstutorials/gromacstutorials-inputs/main/level1/bulk-solution/ff/h2o.itp" target="_blank">h2o.itp</a>

.. |na.itp| raw:: html

    <a href="https://raw.githubusercontent.com/gromacstutorials/gromacstutorials-inputs/main/level1/bulk-solution/ff/na.itp" target="_blank">na.itp</a>

.. |so4.itp| raw:: html

    <a href="https://raw.githubusercontent.com/gromacstutorials/gromacstutorials-inputs/main/level1/bulk-solution/ff/so4.itp" target="_blank">so4.itp</a>

..  container:: justify

    For instance, the *forcefield.itp* file contains a line that specifies
    the combination rules as *comb-rule 2*, which corresponds to
    the well-known Lorentz-Berthelot rule,
    where :math:`\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}` and
    :math:`\sigma_{ij} = (\sigma_{ii}+\sigma_{jj})/2`
    :cite:`lorentz1881ueber,berthelot1898melange`:

..  code-block:: bw

    [ defaults ]
    ; nbfunc  comb-rule  gen-pairs  fudgeLJ  fudgeQQ
    1       2          no         1.0      0.833
    
..  container:: justify

    The fudge parameters specify how the pair interaction between
    fourth neighbors in a residue are handled, which is not relevant for
    the small residues considered here. The *forcefield.itp* file also
    contains the list of atoms, and their respective charge in the units of
    the elementary charge :math:`e`, as well as their respective Lennard-Jones
    parameters :math:`\sigma` (in nanometer)
    and :math:`\epsilon` 
    (in kJ/mol):

..  code-block:: bw

    [ atomtypes ]
    ; name  at.num  mass      charge  ptype  sigma    epsilon
    Na    11      22.9900   1.0000  A      0.23100  0.45000
    OS     8      15.9994  -1.0000  A      0.38600  0.12
    SO    16      32.0600   2.0000  A      0.35500  1.0465
    HW     1       1.0079   0.5270  A      0.00000  0.00000
    OW     8      15.9994   0.0000  A      0.31650  0.77323
    MW     0       0.0000  -1.0540  D      0.00000  0.00000

..  container:: justify

    Here the *ptype* is used to differential the real atoms (A), such as 
    hydrogens and oxygens, from the virtual and massless site of the
    four-point water model (D).

..  container:: justify

    Finally, the *h2o.itp*, *na.itp*, and *so4.itp* files contain information 
    about the residues, such as their exact compositions, which pairs of
    atoms are connected by bonds as well as the parameters for these bonds.
    In the case of the SO\ :sub:`4`\ :sup:`2-`,
    for instance, the sulfur atom forms a bond of equilibrium distance
    :math:`0.152~\text{nm}`
    and rigidity constant :math:`3.7656 \mathrm{e}4 ~ \text{kJ/mol/nm}^2`
    with each of the four oxygen atoms:

..  code-block:: bw

    [ bonds ]
    ;  ai   aj  funct   c0         c1
        1    5    1   0.1520   3.7656e4
        2    5    1   0.1520   3.7656e4
        3    5    1   0.1520   3.7656e4
        4    5    1   0.1520   3.7656e4

3) The input file (.mdp)
------------------------

..  container:: justify

    The input file contains instructions about the
    simulation, such as

    - the number of steps to perform,
    - the thermostat to be used (e.g. Langevin :cite:`schneider1978molecular`, Berendsen :cite:`berendsen1984molecular`),
    - the cut-off for the interactions,
    - the integrator or algorithms (e.g. steepest-descent :cite:`debye1909naherungsformeln`, leap-frog :cite:`allen2017computer`).

..  container:: justify

    In this tutorial, 4 different input files will be
    written to perform respectively an energy
    minimization of the salt solution, an equilibration
    in the NVT ensemble (i.e. with fixed box size), an equilibration
    in the NPT ensemble (i.e. with changing box size), and finally
    a production run.

..  container:: justify

    Input files will be placed in an *inputs/* folder
    that must be created next to *ff/*. 

.. figure:: ../figures/level1/bulk-solution/gromacs_inputs-light.png
    :alt: Gromacs files and structure folder
    :height: 200
    :class: only-light

.. figure:: ../figures/level1/bulk-solution/gromacs_inputs-dark.png
    :alt: Gromacs files and structure folder
    :height: 200
    :class: only-dark

.. container:: figurelegend

    Figure: Structure of the folder with the *.itp*, *.gro*, and *.top*
    files.

..  container:: justify

    The rest of the tutorial focuses on writing the input files and performing the
    molecular dynamics simulation. 

Energy minimization
===================

..  container:: justify

    It is clear from the configuration (.gro) file that the molecules and ions
    are currently in a quite unphysical configuration. It would be risky to
    directly start the molecular dynamics simulation as
    atoms would undergo huge forces and accelerate, and the system could
    eventually explode.

..  container:: justify

    To bring the system into a more favorable state,
    let us perform an energy minimization which
    consists in moving the atoms until the forces between them are reasonable.

..  container:: justify

    Open a blank file, call it *min.mdp*, and save it in the
    *inputs/* folder. Copy the following lines into *min.mdp*:

..  code-block:: bw

    integrator = steep
    nsteps = 5000

..  container:: justify

    These two commands specify to GROMACS that the algorithm
    to be used is the |speepest-descent|
    which moves the atoms following the direction of the largest forces
    until one of the stopping criteria is reached :cite:`debye1909naherungsformeln`. 
    The *nsteps* command specifies the maximum number of steps to perform.

.. |speepest-descent| raw:: html

    <a href="https://manual.gromacs.org/current/reference-manual/algorithms/energy-minimization.html" target="_blank">steepest-descent</a>

..  container:: justify

    To visualize the trajectory of the atoms during the minimization,
    let us also add the following command to the input
    file:

..  code-block:: bw

    nstxout = 10

..  container:: justify

    The *nstxout* parameter requests GROMACS to print the atom
    positions every 10 steps in a *.trr* trajectory file that can
    be read by VMD. We now have a very minimalist input script, let us try
    it. From the terminal, type:

..  code-block:: bash

    gmx grompp -f inputs/min.mdp -c conf.gro -p topol.top -o min -pp min -po min
    gmx mdrun -v -deffnm min

..  container:: justify

    The *grompp* command is used to preprocess the
    files in order to prepare the simulation. The *grompp* command also
    checks the validity of the files. By using
    the *-f*, *-c*, and *-p* keywords, we specify which
    input, configuration, and topology files must be
    used, respectively. The other keywords *-o*, *-pp*, and *-po* are
    used to specify the names of the output that will be
    produced during the run. 
    
..  container:: justify

    The *mdrun* command calls the engine
    performing the computation from the preprocessed
    files (which is recognized thanks to the *-deffnm* keyword). The *-v*
    option enables *verbose* for more
    information printed in the terminal.

..  container:: justify

    If everything works, you should see something like:

..  code-block:: bw

    :-) GROMACS - gmx mdrun, 2023.2 (-:
    Executable:   /usr/bin/gmx
    Data prefix:  /usr
    (...)
    Steepest Descents converged to machine precision in 961 steps,
    but did not reach the requested Fmax < 10.
    Potential Energy  = -4.4659062e+04
    Maximum force     =  2.2164763e+02 on atom 15
    Norm of force     =  4.5976703e+01


..  container:: justify

    The information printed in the terminal indicates us
    that energy minimization has been performed, even
    though the precision that was asked from the default
    parameters were not reached. We can ignore this message,
    as long as the final energy is large and negative,
    the simulation will work just fine. 
    
..  container:: justify

    The final potential energy is large and
    negative, and the maximum force is as small as :math:`220 ~ \text{kJ/mol/nm}`
    (about 0.4 pN). Everything seems alright.
    Let us visualize the trajectories of the atoms during the
    minimization step using VMD by typing in the terminal:

..  code-block:: bash

    vmd conf.gro min.trr

.. figure:: ../figures/level1/bulk-solution/solution-light.webp
    :alt: Gromacs tutorial : Movie showing the motion of the atoms during the energy minimization.
    :class: only-light
    :height: 330

.. figure:: ../figures/level1/bulk-solution/solution-dark.webp
    :alt: Gromacs tutorial : Movie showing the motion of the atoms during the energy minimization.
    :class: only-dark
    :height: 330

.. container:: figurelegend

    Figure: Movie showing the motion of the atoms during the energy minimization.

.. admonition:: Note for VMD users
    :class: info

    You can avoid having
    molecules *cut in half* by the periodic boundary
    conditions by rewriting the trajectory using:

    ..  code-block:: bash
        
        gmx trjconv -f min.trr -s min.tpr -o min_whole.trr -pbc whole


    Then, select the group of your choice for the output.

..  container:: justify

    One can see that the molecules reorient themselves
    into more energetically favorable positions, and that
    the distances between the atoms are being
    progressively homogenized.
    
..  container:: justify

    Let us have a look at the
    evolution of the potential energy of the system. To do
    so, we can use the *gmx energy* command of
    GROMACS. In the terminal, type:

..  code-block:: bash

    gmx energy -f min.edr -o potential-energy-minimization.xvg

..  container:: justify

    Choose *potential* (in my case I have to type *5*),
    then press *Enter* twice. 
    
..  container:: justify

    Here, the portable energy file *min.edr* produced
    by GROMACS during the *minimization run* is used, and the
    result is saved in the *potential-energy-minimization.xvg* file.

.. figure:: ../figures/level1/bulk-solution/potential-energy-min-light.png
    :alt: Gromacs tutorial : plot of the energy versus time.
    :class: only-light

.. figure:: ../figures/level1/bulk-solution/potential-energy-min-dark.png
    :alt: Gromacs tutorial : plot of the energy versus time.
    :class: only-dark

.. container:: figurelegend

    Figure: Evolution of the potential energy :math:`E_\text{p}` as a function of the
    number of steps during energy minimization.

..  container:: justify

    One can see from the energy plot that
    the potential energy is initially huge and positive,
    which is the consequence of atoms being too close from
    one another, as well as molecules being wrongly oriented.
    As the minimization progresses, the potential energy
    rapidly decreases and reaches a large and negative
    value, which is usually the sign
    that the atoms are located at appropriate distances from each other. 

..  container:: justify

    Thank to the energy minimization, the system is now in a favorable state
    and the molecular dynamics simulation can be started safely.

Minimalist NVT input file
=========================

..  container:: justify

    Let us first perform a short (20 picoseconds)
    equilibration in the NVT ensemble. In the NVT ensemble, the number of
    atoms (N) and the volume (V) are maintained fixed, and the
    temperature of the system (T) is adjusted using a thermostat.
    
..  container:: justify

    Let us write a new input script called *nvt.mdp*, and save it in
    the *inputs/* folder. Copy the following lines into it:

..  code-block:: bw

    integrator = md
    nsteps = 20000
    dt = 0.001

..  container:: justify

    Here, the molecular dynamics (md) integrator is used, this is a leap-frog
    algorithm integrating Newton equations of motion. A number of 20000 steps with
    a timestep *dt* equal of :math:`0.001 ~ \text{ps}` will be performed.

..  container:: justify

    Let us also ask GROMACS to print the trajectory in a compressed *xtc* file
    every 1000 steps, or every 1 ps, by adding the following line to
    *nvt.mdp*:

..  code-block:: bw

    nstxout-compressed = 1000

..  container:: justify

    Let us also control the temperature throughout the
    simulation using the v-rescale thermostat, which is
    the Berendsen thermostat with an additional stochastic
    term :cite:`bussi2007canonical`.

..  code-block:: bw

    tcoupl = v-rescale
    ref-t = 360
    tc-grps = system
    tau-t = 0.5

..  container:: justify

    The v-rescale thermostat is known to give
    a proper canonical ensemble. Here, we also specified that the thermostat is
    applied to the entire system using the *tc-grps* option and that the
    damping constant for the thermostat, *tau-t*, is equal to 0.5 ps.

..  container:: justify

    Note that the relatively high temperature of 360 K
    has been chosen here to reduce the viscosity of the solution and
    decrease the equilibration duration. 
    
..  container:: justify
    
    We now have a minimalist input file for performing
    the first NVT simulation. Run it by typing in the terminal:

..  code-block:: bw

    gmx grompp -f inputs/nvt.mdp -c min.gro -p topol.top -o nvt -pp nvt -po nvt
    gmx mdrun -v -deffnm nvt

..  container:: justify

    Here *-c min.gro* ensures that the previously
    minimized configuration is used as a starting point.
    
..  container:: justify

    After the completion of the simulation, we can
    ensure that the system temperature indeed reached
    the value of 360 K by using the energy command of
    GROMACS. In the terminal, type:

..  code-block:: bw

    gmx energy -f nvt.edr -o temperature-nvt-minimal.xvg

..  container:: justify

    and choose *temperature*.
    
..  container:: justify

    From the generated *temperature-nvt-minimal.xvg* file, one can see that temperature
    started from 0 K, which was expected since the atoms have no velocity
    during a minimization step, and reaches a temperature slightly larger than the
    requested 360 K after a duration of a few picoseconds.

..  container:: justify

    In general, it is better to perform a longer equilibration, but simulation
    durations are kept as short as possible for these tutorials. 

.. figure:: ../figures/level1/bulk-solution/temperature-nvt-minimal-light.png
    :alt: Gromacs tutorial : temperature versus time.
    :class: only-light

.. figure:: ../figures/level1/bulk-solution/temperature-nvt-minimal-dark.png
    :alt: Gromacs tutorial : temperature versus time.
    :class: only-dark

.. container:: figurelegend

    Figure: Evolution of the temperature :math:`T` as a function of
    the time :math:`t` during the NVT equilibration. The dashed line is the
    requested temperature of 360 K.

Improving the NVT input
=======================

..  container:: justify

    So far, very few commands have been placed in the
    *.mdp* input file, meaning that most of the instructions
    have been taken by GROMACS from the default
    parameters. You can find what parameters were used
    during the last nvt run by opening the new *nvt.mdp*
    file that has been created in the main folder.
    Exploring this new *nvt.mdp* file shows us that, for
    instance, plain cut-off Coulomb interactions have
    been used:

..  code-block:: bw

    (...)
    ; Method for doing electrostatics
    coulombtype = Cut-off
    (...)

..  container:: justify

    For this system, computing the long-range Coulomb interactions is necessary,
    because electrostatic forces between charged particles decay slowly,
    as the inverse of the square of the distance between them, :math:`1/r^2`.
    
..  container:: justify

    In addition, the thermostating of the system should be improved, given that
    the temperature of the system is slightly larger than the desired temperature.
    For instance, separate thermostats can be applied to the water molecules
    and the ions.

..  container:: justify

    Let us improve the input used for the NVT step.
    First, in the *nvt.mdp* file, let us impose the calculation of long-range
    electrostatic, by the use of the
    long-range fast smooth particle-mesh ewald (SPME)
    electrostatics with Fourier spacing of :math:`0.1~\text{nm}`, order
    of 4, and cut-off of :math:`1~\text{nm}`
    :cite:`darden1993particle, essmann1995smooth`:

..  code-block:: bw

    coulombtype = pme
    fourierspacing = 0.1
    pme-order = 4
    rcoulomb = 1.0

..  container:: justify

    Here, the cut-off *rcoulomb* separates the short-range interactions from the
    long-range interactions. Long-range interactions are treated in the
    reciprocal space, while short-range interactions are computed directly.

..  container:: justify

    Let us also impose how the short-range van der Waals interactions
    should be treated by GROMACS, as well as the cut-off *rvdw*
    of :math:`1~\text{nm}`:

..  code-block:: bw

    vdw-type = Cut-off
    rvdw = 1.0

..  container:: justify

    Let us use the LINCS algorithm to constrain the hydrogen
    bonds :cite:`hess1997lincs`. The water
    molecules will thus be treated as rigid, which is generally better given
    the fast vibration of the hydrogen bonds.

..  code-block:: bw

    constraint-algorithm = lincs
    constraints = hbonds
    continuation = no

..  container:: justify

    Let us also use separate temperature baths for
    the water molecules and the ions. Here, the ions are included
    in the default GROMACS group called *non-water*.
    Within *nvt.mdp*, replace the following lines:

..  code-block:: bw

    tcoupl = v-rescale
    ref-t = 360
    tc-grps = system
    tau-t = 0.5

..  container:: justify

    by:

..  code-block:: bw

    tcoupl = v-rescale
    tc-grps = Water non-Water
    tau-t = 0.5 0.5
    ref-t = 360 360

..  container:: justify

    Now, the same temperature :math:`T = 360 ~ \text{K}` is imposed to the
    two groups with the same characteristic time :math:`\tau = 0.5 ~ \text{ps}`.
    
..  container:: justify
  
    Let us also specify the neighbor searching parameters:

..  code-block:: bw

    cutoff-scheme = Verlet
    nstlist = 10
    ns_type = grid

..  container:: justify

    Let us give an initial kick to the atom so that the initial
    total velocities give the desired temperature of 360 K instead of 0 K
    as previously:

..  code-block:: bw

    gen-vel = yes
    gen-temp = 360

..  container:: justify

    Finally, let us cancel the translational of the center of mass
    of the system:

..  code-block:: bw

    comm_mode = linear
    comm_grps = system

..  container:: justify

    Run again GROMACS using this new input script. One difference with the
    previous (minimalist) NVT run is
    the temperature at the beginning of the run. The final
    temperature is also much closer to the desired temperature of 360 K.

.. figure:: ../figures/level1/bulk-solution/temperature-nvt-light.png
    :alt: Gromacs tutorial : temperature versus time.
    :class: only-light

.. figure:: ../figures/level1/bulk-solution/temperature-nvt-dark.png
    :alt: Gromacs tutorial : temperature versus time.
    :class: only-dark

.. container:: figurelegend

    Figure: Evolution of the temperature as a function of the time
    during the NVT equilibration.

Adjust the density using NPT
============================

..  container:: justify

    Now that the temperature of the system is properly equilibrated,
    let us continue the simulation using the
    NPT ensemble, where the pressure of the system is imposed by a barostat
    and the volume of the box is allowed to relax. During NPT relaxation, the
    density of the fluid should converge toward its equilibrium value.
    Create a new input script, call it *npt.mdp*, and
    copy the following lines in it:

..  code-block:: bw

    integrator = md
    nsteps = 80000
    dt = 0.001

    comm_mode = linear
    comm_grps = system

    cutoff-scheme = Verlet
    nstlist = 10
    ns_type = grid

    nstlog = 100
    nstenergy = 100
    nstxout-compressed = 1000

    vdw-type = Cut-off
    rvdw = 1.0

    coulombtype = pme
    fourierspacing = 0.1
    pme-order = 4
    rcoulomb = 1.0

    constraint-algorithm = lincs
    constraints = hbonds

    tcoupl = v-rescale
    ld-seed = 48456
    tc-grps = Water non-Water
    tau-t = 0.5 0.5
    ref-t = 360 360

    pcoupl = C-rescale
    Pcoupltype = isotropic
    tau_p = 1.0
    ref_p = 1.0
    compressibility = 4.5e-5

..  container:: justify

    The main difference with the previous NVT script is
    the addition of the isotropic C-rescale pressure
    coupling with a target pressure of 1 bar :cite:`bernetti2020pressure`.
    Another difference is the addition of the *nstlog* and
    *nstenergy* commands to control the frequency at
    which information is printed in the log file and in
    the energy file (*edr*). Note also the removing the
    *gen-vel* commands, because the atoms already have a velocity. 
    
..  container:: justify

    Run the NPT equilibration using:

..  code-block:: bash 

    gmx grompp -f inputs/npt.mdp -c nvt.gro -p topol.top -o npt -pp npt -po npt
    gmx mdrun -v -deffnm npt

..  container:: justify

    Let us have a look a the temperature, the pressure, and the
    volume of the box during the NPT step using the *gmx energy*
    command 3 consecutive times:

..  code-block:: bash

    gmx energy -f npt.edr -o temperature-npt.xvg
    gmx energy -f npt.edr -o pressure-npt.xvg
    gmx energy -f npt.edr -o density-npt.xvg

..  container:: justify

    Choose respectively *temperature*, *pressure* and *density*.
    This is what I see:

.. figure:: ../figures/level1/bulk-solution/temperature-npt-light.png
    :alt: Gromacs tutorial : NPT equilibration
    :class: only-light

.. figure:: ../figures/level1/bulk-solution/temperature-npt-dark.png
    :alt: Gromacs tutorial : NPT equilibration
    :class: only-dark

..  container:: figurelegend

    Figure: Evolution of the temperature :math:`T` (a),
    pressure :math:`p` (b),
    and fluid density :math:`\rho` (c) as a
    function of the time during the NPT equilibration.

..  container:: justify

    The results show that the temperature remains well
    controlled during the NPT run, and
    that the fluid density was initially too small,
    i.e. :math:`\rho \approx 600\,\mathrm{kg}/\mathrm{m}^3`.
    Due to the change in volume induced by the barostat, the fluid density 
    gently reaches its equilibrium value of about :math:`1000\,\mathrm{kg}/\mathrm{m}^3`
    after approximately 40 pico-seconds. Once the system has reached its equilibrium
    density, the pressure stabilizes itself near the desired value of 1 bar.
    
..  container:: justify

    The pressure curve reveals large oscillations in the pressure, with the pressure alternating between large negative
    values and large positive values. These large oscillations are typical in molecular dynamics, and not a source of
    concern here.

Radial distribution function
============================

..  container:: justify

    Let us perform a :math:`400~\text{ps}` run in the NVT ensemble, during which the atom positions will be printed every
    pico-second. The trajectory will then be used to measure radial distribution functions and probe the solvation
    environment of the ions. 
    
..  container:: justify
    
    Create a new input file within the *inputs/* folder, call it *production.mdp*, and copy the following lines into it:

..  code-block:: bw 

    integrator = md
    nsteps = 400000
    dt = 0.001

    comm_mode = linear
    comm_grps = system

    cutoff-scheme = Verlet
    nstlist = 10
    ns_type = grid

    nstlog = 100
    nstenergy = 100
    nstxout-compressed = 1000

    vdw-type = Cut-off
    rvdw = 1.0

    coulombtype = pme
    fourierspacing = 0.1
    pme-order = 4
    rcoulomb = 1.0

    constraint-algorithm = lincs
    constraints = hbonds

    tcoupl = v-rescale
    ld-seed = 48456
    tc-grps = Water non-Water
    tau-t = 0.5 0.5
    ref_t = 360 360

..  container:: justify

    Run it using:

..  code-block:: bash 

    gmx grompp -f inputs/production.mdp -c npt.gro -p topol.top -o production -pp production -po production
    gmx mdrun -v -deffnm production

..  container:: justify

    When the simulation is completed, let us compute the radial
    distribution functions between :math:`\text{Na}^+` and
    :math:`\text{H}_2\text{O}`, 
    :math:`\text{SO}_4^{2-}` and 
    :math:`\text{H}_2\text{O}`,
    as well as in between :math:`\text{H}_2\text{O}` molecules.
    This can be done using the *gmx rdf* command as follows:
    
..  code-block:: bash 

    gmx rdf -f production.xtc -s production.tpr -o na-sol-rdf.xvg

..  container:: justify

    Selecting the sodium ions, and then the water. Repeat the same operation for 
    the sulfate and water, and for the water and water. For the water-water
    RDF, it is better to exclude the intra-molecular contribution using
    the *-excl* option, as follows:

..  code-block:: bash 

    gmx rdf -f production.xtc -s production.tpr -o sol-sol-rdf.xvg -excl
    
.. figure:: ../figures/level1/bulk-solution/rdf-production-light.png
    :alt: Gromacs tutorial RDF radial distribution function
    :class: only-light

.. figure:: ../figures/level1/bulk-solution/rdf-production-dark.png
    :alt: Gromacs tutorial RDF radial distribution function
    :class: only-dark

.. container:: figurelegend

    Figure: Radial distribution functions (RDF) as calculated between sodium
    and water, between sulfate and water, and finally between water and water.
    
..  container:: justify

    The radial distribution functions highlight the typical distance between
    the different species of the fluid. For instance, it can be seen that
    there is a strong hydration layer of water around sodium ions at a typical
    distance of :math:`2.4 ~ \text{Å}` from the center of the sodium ion.

Mean square displacement
========================

..  container:: justify

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

..  code-block:: bash 

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

..  container:: justify

    Select the :math:`\text{SO}_4^{2-}` ions (in my case, it is done by typing
    *2*), and then press *ctrl D*. Repeat the same operation for the sodium
    ions and for the water molecules.
    
..  container:: justify

    The slope of the MSD in the limit of long times gives an estimate of the diffusion
    coefficient, following :math:`D = \text{MSD} / 2 d t`,
    where :math:`d = 3` is the dimension of the system. Here,
    I find a value of :math:`1.4 \mathrm{e}-5 ~ \text{cm}^2/\text{s}` for the
    diffusion coefficient of the sulfur ions. 

..  container:: justify
    
    Repeat the same for :math:`\text{Na}^+` and water. 
    
..  container:: justify

    For sodium, I find a value of :math:`1.6 \mathrm{e}-5 ~ \text{cm}^2/\text{s}`
    for the diffusion coefficient,
    and for water :math:`5.3 \mathrm{e}-5 ~ \text{cm}^2/\text{s}`.
    For comparison, the experimental diffusion coefficient of *pure* water at 
    temperature :math:`T = 360~\text{K}`
    is :math:`7.3 \mathrm{e}-5 ~ \text{cm}^2/\text{s}` :cite:`simpson1958diffusion`.
    In the presence of ions, the diffusion coefficient of water is expected to
    be reduced.

.. figure:: ../figures/level1/bulk-solution/msd-production-light.png
    :alt: Gromacs tutorial : diffusion coefficient
    :class: only-light

.. figure:: ../figures/level1/bulk-solution/msd-production-dark.png
    :alt: Gromacs tutorial : diffusion coefficient
    :class: only-dark

.. container:: figurelegend

    Figure: Mean square displacement (msd) for the three species. The dashed line
    highlight the proportionality between msd and time :math:`t` which is expected
    at long times, when the system reaches the diffusive regime.

.. admonition:: About diffusion coefficient measurement in molecular simulations
    :class: info

    In principle, diffusion coefficients obtained from the MSD in a
    finite-sized box must be corrected, but this is beyond the scope of
    the present tutorial :cite:`loche2021transferable`.

.. include:: ../../non-tutorials/accessfile.rst

.. 
    Going further with exercises
    ============================

    ..  container:: justify

        The solutions can be found in the GitHub repository of GROMACS tutorials.

    Calculate the structure factor
    ------------------------------

    ..  container:: justify

        Use the *gmx scattering* command and extract its structure factor, or :math:`S(q)`.
        The structure factor is a fundamental quantity that it commonly measured
        experimentally through Small Angle X-ray Scattering (SAXS) or Small-Angle
        Neutron Scattering (SANS) :cite:`hansen2013theory, sedlmeier2011spatial`.