.. _solvation-energy-label: Solvation energy **************** .. container:: hatnote Calculating the free energy of solvation of a graphene-like molecule .. figure:: ../figures/level3/solvation-energy/video-HBC-light-2.webp :alt: HBC (graphene-like) molecule in water :class: only-light :height: 250 :align: right .. figure:: ../figures/level3/solvation-energy/video-HBC-dark-2.webp :alt: HBC (graphene-like) molecule in water :class: only-dark :height: 250 :align: right The objective of this tutorial is to use GROMACS to perform a molecular simulation of a large molecule in water. By progressively switching off the interactions between the molecule and water, the free energy of solvation will be calculated. The large and flat molecule used here is a graphene-like and discoid molecule named hexabenzocoronene (HBC) of formula :math:`\text{C}_{42}\text{H}_{18}`. The TIP4P/epsilon model is used for the water :cite:`fuentes-azcatlNonPolarizableForceField2014`. .. include:: ../../non-tutorials/recommand-salt.rst .. include:: ../../non-tutorials/needhelp.rst .. include:: ../../non-tutorials/GROMACS2024.2.rst Input files =========== Create two folders side-by side. Name the two folders *preparation/* and *solvation/*, and go to *preparation/*. Download the configuration files for the HBC molecule by clicking |FJEW_allatom_optimised_geometry.pdb|, and save it in the *preparation/* folder. The molecule was downloaded from the |atb-HBC| on the automated topology builder (ATB) :cite:`maldeAutomatedForceField2011`. .. |FJEW_allatom_optimised_geometry.pdb| raw:: html this page .. |atb-HBC| raw:: html atb Create the configuration file ----------------------------- First, let us convert the *.pdb* file into a *.gro* file within a cubic box of lateral size 3 nanometers using the *gmx trjconv* command. Type the following command in a terminal: .. code-block:: bw gmx trjconv -f FJEW_allatom_optimised_geometry.pdb -s FJEW_allatom_optimised_geometry.pdb -o hbc.gro -box 3 3 3 -center Select *system* for both centering and output. .. figure:: ../figures/level3/solvation-energy/hbc-light.png :alt: Gromacs initial configuration of HBC graphene molecule :class: only-light :height: 250 :align: center .. figure:: ../figures/level3/solvation-energy/hbc-dark.png :alt: Gromacs initial configuration of HBC graphene molecule :class: only-dark :height: 250 :align: center .. container:: figurelegend Figure: HBC molecule as seen with VMD with carbon atoms in gray and hydrogen atoms in white. The honeycomb structure of the HBC is similar to that of graphene. Alternatively, you can download the |solvation-hbc.gro| I have generated, and continue with the tutorial. .. |solvation-hbc.gro| raw:: html hbc.gro Create the topology file ------------------------ Create a folder named *ff/* within the *preparation/* folder. Copy the force field parameters from the following *zip* file by clicking |ff-itp.zip|. Both *FJEW_GROMACS_G54A7FF_allatom.itp* file and *gromos54a7_atb.ff/* folder were downloaded from the |atb-HBC|. .. |ff-itp.zip| raw:: html here Then, let us write the topology (*top*) file by simply creating a blank file named *topol.top* within the *preparation/* folder. Copy the following lines into *topol.top*: .. code-block:: bw #include "ff/gromos54a7_atb.ff/forcefield.itp" #include "ff/FJEW_GROMACS_G54A7FF_allatom.itp" [ system ] Single HBC molecule [ molecules ] FJEW 1 Solvate the HBC in water ------------------------ Let us add water molecules to the system. First, download the tip4p water configuration (*.gro*) file |tip4p.gro|, and copy it in the *preparation/* folder. This is a fourth point water model with additional massless site where the charge of the oxygen atom is placed. Then, in order to add tip4p water molecules to both *.gro* and *.top* file, use the *gmx solvate* command as follows: .. |tip4p.gro| raw:: html here .. code-block:: bw gmx solvate -cs tip4p.gro -cp hbc.gro -o solvated.gro -p topol.top The new *solvated.gro* file contains all *8804* atoms from the HBC molecule (called FJEW) and the water molecules. A new line *SOL 887* also appeared in the topology *.top* file: .. code-block:: bw [ molecules ] FJEW 1 SOL 887 Alternatively, you can download the *solvated.gro* file I have generated by clicking |solvated.gro|, and continue with the tutorial. .. |solvated.gro| raw:: html here Finally, save the topology file for the water, the |h2o.itp| file, in the *ff/* folder and add the *#include "ff/h2o.itp"* line to the *topol.top* file: .. |h2o.itp| raw:: html h2o.itp .. code-block:: bw #include "ff/gromos54a7_atb.ff/forcefield.itp" #include "ff/FJEW_GROMACS_G54A7FF_allatom.itp" #include "ff/h2o.itp" System equilibration ==================== The system is now ready for the simulations. Let us first equilibrate it before measuring the solvation energy of the HBC molecule. Energy minimization ------------------- Create an *inputs/* folder inside the *preparation/* folder, and create a new blank file called *min.mdp* in it. Copy the following lines into *min.mdp*: .. code-block:: bw integrator = steep nsteps = 50 nstxout = 10 cutoff-scheme = Verlet nstlist = 10 ns_type = grid couple-intramol=yes vdw-type = Cut-off rvdw = 1.2 coulombtype = pme fourierspacing = 0.1 pme-order = 4 rcoulomb = 1.2 constraint-algorithm = lincs constraints = hbonds All these lines have been seen in the previous tutorials. In short, with this input script, GROMACS will perform a steepest descent by updating the atom positions according to the largest forces directions, until the energy and maximum force reach a reasonable value. Apply the minimization to the solvated box using *gmx grompp* and *gmx mdrun*: .. code-block:: bash gmx grompp -f inputs/min.mdp -c solvated.gro -p topol.top -o min -pp min -po min -maxwarn 1 gmx mdrun -v -deffnm min Here, the *-maxwarn 1* option is used to ignore a WARNING from GROMACS about some issue with the force field. For this tutorial, we can safely ignore this WARNING. Let us visualize the atom trajectories during the minimization step using VMD by typing: .. code-block:: bash vmd solvate.gro min.trr .. figure:: ../figures/level3/solvation-energy/minimize-light.png :alt: Gromacs HBC (graphene) molecule after minimization in water :class: only-light :height: 400 :align: center .. figure:: ../figures/level3/solvation-energy/minimize-dark.png :alt: Gromacs HBC (graphene) molecule after minimization in water :class: only-dark :height: 400 :align: center .. container:: figurelegend Figure: The system after energy minimization showing the HBC molecule in water. The water is represented as a transparent field. NVT and NPT equilibration ------------------------- Let us perform successively a NVT and a NPT relaxation steps. Copy the |solvation-nvt.mdp| and the |solvation-npt.mdp| files into the inputs folder, and run them both using: .. |solvation-nvt.mdp| raw:: html nvt.mdp .. |solvation-npt.mdp| raw:: html npt.mdp .. code-block:: bash gmx grompp -f inputs/nvt.mdp -c min.gro -p topol.top -o nvt -pp nvt -po nvt -maxwarn 1 gmx mdrun -v -deffnm nvt gmx grompp -f inputs/npt.mdp -c nvt.gro -p topol.top -o npt -pp npt -po npt -maxwarn 1 gmx mdrun -v -deffnm npt .. figure:: ../figures/level3/solvation-energy/nvtnpt-light.webp :alt: Gromacs HBC (graphene) molecule during NVT and NPT equilibration :class: only-light :height: 400 :align: center .. figure:: ../figures/level3/solvation-energy/nvtnpt-dark.webp :alt: Gromacs HBC (graphene) molecule during NVT and NPT equilibration :class: only-dark :height: 400 :align: center .. container:: figurelegend Figure: Movie showing the motion of the atoms during the NVT and NPT equilibration steps. For clarity, the water molecules are represented as a continuum field. Solvation energy measurement ============================ Files preparation ----------------- The equilibration of the system is complete. Let us perform the solvation free energy calculation, for which 21 independent simulations will be performed. .. admonition:: About free energy calculation :class: info The interactions between the HBC molecule and the water are progressively turned-off, thus effectively mimicking the HBC molecule moving from bulk water to vacuum. Then, the free energy difference between the fully solvated and the fully decoupled configurations is measured. Within the *solvation/* folder, create an *inputs/* folders, and copy the following input file into it: |equilibration.mdp|. .. |equilibration.mdp| raw:: html equilibration.mdp This input file starts similarly as the inputs previously used in this tutorial, with the exception of the integrator. Instead of the *md* integrator, a *sd* for stochastic dynamics integrator: .. code-block:: bw integrator = sd nsteps = 20000 dt = 0.001 This stochastic integrator creates a Langevin dynamics by adding a friction and a noise term to Newton equations of motion. The *sd* integrator also serves as a thermostat, therefore *tcoupl* is set to *No*: .. code-block:: bw tcoupl = No ld-seed = 48456 tc-grps = Water non-Water tau-t = 0.5 0.5 ref-t = 300 300 A stochastic integrator can be a better option for free energy measurements, as it generates a better sampling, while also imposing a strong control of the temperature. This is particularly useful when the molecules are completely decoupled. The rest of the input deals with the progressive decoupling of the HBC molecule from the water: .. code-block:: bw free_energy = yes vdw-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 coul-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 The *vdw-lambdas* is used to turn-off (when *vdw-lambdas = 0*) or turn-on (when *vdw-lambdas = 1*) the van der Waals interactions, and the *coul-lambdas* is used to turn-off/turn-on the Coulomb interactions. Here, there are 21 possible states, from *vdw-lambdas = coul-lambdas = 0.0*, where the HBC molecule is fully decoupled from the water, to *vdw-lambdas = coul-lambdas = 1.0*, where the HBC molecule is fully coupled to the water. All the other state correspond to situation where the HBC molecule is partially coupled with the water. The *sc-alpha* and *sc-power* options are used to turn-on a soft core repulsion between the HBC and the water molecules. This is necessary to avoid overlap between the decoupled atoms. The *init-lambda-state* is an integer that specifies which values of *vdw-lambdas* and *coul-lambdas* is used. When *init-lambda-state*, then *vdw-lambdas = coul-lambdas = 0.0*. 21 one independent simulations will be performed, each with a different value of init-lambda-state: .. code-block:: bw init-lambda-state = 0 The *couple-lambda0* and *couple-lambda1* are used to specify that indeed interaction are turn-off when the *lambdas* are 0, and turn-on when the *lambdas* as 1: .. code-block:: bw couple-lambda0 = none couple-lambda1 = vdw-q The parameter *nstdhdl* controls the frequency at which information are being printed in a *.xvg* file during the production run. .. code-block:: bw nstdhdl = 0 For the equilibration, there is no need of printing information. For the production runs, a value :math:`>0` will be used for nstdhdl. The 2 last lines impose that lambda points will be written out, and which molecule will be used for calculating solvation free energies: .. code-block:: bw calc_lambda_neighbors = -1 couple-moltype = FJEW Let us create a second input file almost identical to *equilibration.mdp*. Duplicate *equilibration.mdp*, name the duplicated file *production.mdp*. Within *production.mdp*, simply change *nstdhdl* from 0 to 100, so that information about the state of the system will be printed by GROMACS every 100 step during the production runs. Finally, copy the previous *topol.top* file from the *preparation/* folder into *solvation/* folder. In oder to avoid duplicating the force field folder *ff/*, modify the path to the *.itp* files as follow: .. code-block:: bw #include "../preparation/ff/gromos54a7_atb.ff/forcefield.itp" #include "../preparation/ff/FJEW_GROMACS_G54A7FF_allatom.itp" #include "../preparation/ff/h2o.itp" [ system ] Single HBC molecule in water [ molecules ] FJEW 1 SOL 887 Perform a 21-step loop ---------------------- We need to a total of 2 x 21 simulations, 2 simulations per value of *init-lambda-state*. This can be done using a small bash script with the *sed* (for *stream editor*) command. Create a new bash file fine within the 'solvation/' folder, call it *local-run.sh* can copy the following lines in it: .. code-block:: bw #/bin/bash set -e # create folder for analysis mkdir -p dhdl # loop on the 21 lambda state for state in $(seq 0 20) do # create folder for the lambda state DIRNAME=lambdastate${state} mkdir -p $DIRNAME # update the value of init-lambda-state newline='init-lambda-state = '$state linetoreplace=$(cat inputs/equilibration.mdp | grep init-lambda-state) sed -i '/'"$linetoreplace"'/c\'"$newline" inputs/equilibration.mdp linetoreplace=$(cat inputs/production.mdp | grep init-lambda-state) sed -i '/'"$linetoreplace"'/c\'"$newline" inputs/production.mdp gmx grompp -f inputs/equilibration.mdp -c ../preparation/npt.gro -p topol.top -o equilibration -pp equilibration -po equilibration -maxwarn 1 gmx mdrun -v -deffnm equilibration -nt 4 gmx grompp -f inputs/production.mdp -c equilibration.gro -p topol.top -o production -pp production -po production -maxwarn 1 gmx mdrun -v -deffnm production -nt 4 mv production.* $DIRNAME mv equilibration.* $DIRNAME # create links for the analysis cd dhdl/ ln -sf ../$DIRNAME/production.xvg md$state.xvg cd .. done Within this bash script, the variable *state* increases from 0 to 20 in a for loop. At eash step of the loop, a folder *lambdastateX* is created, where *X* goes from 0 to 20. Then, the *sed* command is used twice to update the value of *init-lambda-state* in both *equilibration.mdp* and *production.mdp*. Then, GROMACS is used to run the *equilibration.mdp* script, and then the *production.mdp* script. When the simulations are done, the generated files are moved into the *lambdastateX* folder. Finally the *ln* command creates a link toward the *production.xvg* file within the *dhdl/* folder. Execute the bash script by typing: .. code-block:: bash bash createfolders.sh Completing the 2 x 21 simulations may take a while, depending on your computer. When the simulation is complete, go the dhdl folder, and call the *gmx bar* command: .. code-block:: bash gmx bar -f *.xvg The value of the solvation energy is printed in the terminal. In my case, I see: .. code-block:: bash DG -64.62 +/- 6.35 This indicate that the solvation energy is of :math:`-64.6 \pm 6.3~\text{kJ/mol}`. Note however that the present simulations are too short to give a reliable result. To accurately measure the solvation energy of a molecule, you should use much longer equilibration, typically one nanosecond, as well as much longer production runs, typically several nanoseconds. .. include:: ../../non-tutorials/accessfile.rst