Tutorial: Using GromPy in GCMC mode - GromPy/GromPy GitHub Wiki

This tutorial will let the user undergo all necessary steps to perform simulations in the grand-canonical ensemble using the GROMACS molecular simulation library via the GromPy interface. After going through this tutorial you will be able to calculate an equation of state (EOS) of the Lennard-Jones (LJ) fluid. The system consists of mono-atomic water molceules defined in the MARTINI force field. The intermolecular interactions between such molecules are modelled using the Lennard-Jones potential only. We assume that you use bash shell of the Linux operating system. If you are unfamiliar with some Linux commands, you can find documentation in the corresponding man pages or on the web.

1. Downloading and installing the required files

  • Go to your favorite working directory: $WORK, e.g. WORK=/home/user/Simulation
    • cd $WORK
  • Use git to download the GromPy files
    • git clone git://github.com/GromPy/GromPy.git
  • Go to the GromPy directory
    • cd GromPy
  • Compile the GROMACS source code
    • cd gromacs
    • cd gromacs-4.0.7
    • tar -xzf gromacs-4.0.7-git.tar.gz
    • cd gromacs-4.0.7-git
    • patch -p1 < ../grompy_4.0.7_patch.diff
    • ./bootstrap
    • ./configure --enable-shared --enable-grompy --prefix=$PWD/install CFLAGS="-O2 -fPIC"
    • make -j 10
    • make install
  • Save the GROMACS install path to a variable
    • GMXINSTALLDIR=$PWD/install
  • Go to the tutorial directory
    • cd ../../../GcmcTutorial
  • We now need to source the GROMACS and GromPy environment variables
    THIS NEEDS TO BE DONE WHENEVER YOU USE GromPy IN A NEWLY OPENED SHELL (SO REMEMBER THE GROMACS INSTALLATION DIRECTORY)!!
    • source $GMXINSTALLDIR/bin/GMXRC
    • source ./SourceGromPyEnv.sh $GMXINSTALLDIR

2. Calculation of the Lennard-Jones equation of state

  • The initial starting structure of the LJ system is an equilibrated one containing N = 400 LJ particles in a cubic box. Since it is equilibrated, there is no need to energy minimize the system. We will however equilibrate it using MD for the desired temperature T = 773 K, just to be sure. This is done in two steps: (1) equilibration using a small time step and (2) production at the target time step dt = 0.02 ps.
  • cd GcmcLj
  • mkdir mdeq
  • cd mdeq
  • grompp -f ../mdp/mdeq.mdp -c ../initstructure/water.gro -p ../top/W400.top
  • mdrun
  • cd ..
  • mkdir mdprod
  • cd mdprod
  • grompp -f ../mdp/mdprod.mdp -c ../mdeq/confout.gro -p ../top/W400.top -maxwarn 1
  • mdrun
  • To enable sampling of LJ particle numbers in the range N ∈ [0, 450], we need to generate 450 .tpr files (for N = 0 we do not need to perform energy calculations ;-)). We delete LJ particles from the equilibrated 400.gro structure to generate the N ∈ [1, 399] range of .gro files and then we add LJ particles for the N ∈ [401, 450] range of .gro files. Remember that GromPy only needs an equilibrated starting structure (one that does not contain particle overlaps). So as long as we start GromPy using a N ∈ [1, 400] structure, it will run fine. For the higher N range we can just generate structures by adding a new LJ particle the last particle on top of e.g. the last particle. Before generating the .tpr files, we also need to generate 450 topology (.top) files.
    • cd ../gro
    • ln -sf ../mdprod/confout.gro ./W400.gro
    • ./GenerateStartingStructures.sh 1 450 W400.gro
    • cd ../top
    • ./GenerateTopologies.sh 1 450 W400.top
    • cd ../tpr
    • ./GenerateTprs.sh 1 450
  • We are now ready to perform grand-canonical simulations. In this tutorial we will generate an equation of state in the ÎŒVT ensemble at T = 773 K
    • cd ../
    • ./RunEos.sh

3. Analysis

  • After performing the simulations, we will analyze the results by generating block averages of the results. We post process the resulting file to obtain the equation of state: the excess chemical potential as a function of density.
    • ./DoBlockAveraging.sh
  • Use your favorite plotting tool to check out the results in file rho of mu with units.dat. If all works well, the results should agree with the NVT equation of state at T = 773 K.