GROMACS Tutorial 6 -- Excess Chemical Potential of Methane using Test Particle Insertion
In this tutorial, we'll be using test particle insertion (TPI) to calculate the excess chemical potential of methane solvation in water. Most users are unaware that GROMACS has a built-in method for running TPI. This tutorial will not be a comprehensive discussion on the statistical mechanics of TPI, but will address issues when needed. The user is encouraged to seek out scientific resources regarding this method.
TPI involves perturbing some states to some other, very similar states. We will be taking bulk water and inserting a methane particle and measuring the potential energy change from this. There is a statistical mechanical relationship between this change in potential energy and the excess chemical potential. For us, state A is the bulk water system, and state B is the water system with a methane.
With GROMACS you need to run state A as a normal MD simulation. We already did this for our case of bulk water in Tutorial 1. We'll reuse the output trajectory files for inserting the methane.
Create water system
Follow Tutorial 1 to run a system containing TIP4PEW water.
Add test particle to topology file
Our original topology file just had water. In the new topology file, we simply
need to add 1 test particle, and it needs to be the last molecule in the system.
opls_066 for the particle's atom type which is OPLS's united atom
methane. Here's what my final topology file looks like (the number of waters
will be different for your system):
#include "oplsaa.ff/forcefield.itp" #include "oplsaa.ff/tip4pew.itp" [ moleculetype ] ; Name nrexcl Methane 3 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 opls_066 1 CH4 C 1 0 16.043 [ System ] Methane in water [ Molecules ] SOL 395 Methane 1
Add test particle to gro file
You also need to add the test particle to the gro file. Simply edit
(or any of the other
.gro files uses) and add a line at the end containing the
test particle's position (right before the box coordinates). The line I added
looks like this:
396CH4 C 1581 0.000 0.000 0.000
The actual position doesn't matter; GROMACSS just wants a placeholder for the
test particle. Additionally, you need to add 1 to the total number of particles
in the system on the second line of the
We only need one parameter file for TPI. Simply copy
prd.mdp from your bulk
water simulation and change
tpi. You should change
the number of insertions per frame that you want to attempt. I chose
steps for my simulation. You will also need
Verlet has not been implemented for
For the simulation, we are just rerunning the bulk water simulation using the
saved trajectory file (which was named
prd.xtc in the first tutorial). To do
this first run
$ gmx grompp -f mdp/tpi.mdp -o tpi.tpr -po tpi.mdp -pp tpi.top -c conf.gro
Now use the
-rerun flag with
$ gmx mdrun mdrun -s tpi.tpr -o tpi.trr -x tpi.xtc -c tpi.gro -e tpi.edr -g tpi.log -rerun prd.xtc
The log file, named
tpi.log in this case, contains a line with the
average volume and the average excess chemical potential. My two lines looked
<V> = 1.18704e+01 nm^3 <mu> = 8.81230e+00 kJ/mol
<mu> is output in kJ/mol, but if we convert it to kcal/mol we get 2.106
kcal/mol. This is in line with our results from the free energy of solvation
done in Tutorial 4
using the lambda-coupling method where I got 2.289 kcal/mol. The difference can
be attributed to the usage of an all-atom model with the free energy of
solvation simulations and a united-atom model in this case.
In this tutorial, we looked at how to use GROMACS to perform test particle insertion in order to get the excess chemical potential of a united-atom OPLS methane.
Author: Wes Barnett, Vedran Miletić