Using energy groups in Gromacs

Motivation

Computing energies of groups within a simulation is useful for measuring energies of a specific molecule or testing changes in forcefield parameters. In this tutorial, we will go through computing Lennard-Jones (LJ) energies of a particular group in a simulation from Gromacs, summarized below:

Goals of this tutorial

At the end of this tutorial, you should be able to:

  • Visualize the system using visual molecular dynamics (VMD) studio

  • Generate energy groups from molecular dynamics parameter (mdp) files

  • Using the energy groups, compute Lennard-Jones energies using "gmx -rerun" command

Software required for this tutorial

Files for this tutorial

All input/output files are downloadable through GitLab [Link]. To download, click the "Download" button and select "Download this directory".

  • P_in_S-PET_trimer-acetone-300_K: Simulation folder containing trajectories.

    • "poly_solv_prod*" prefix contains 10 ns production trajectories in the NPT ensemble

    • poly_solv.top: topology file for the system

    • prod.mdp: mdp file for production run

  • images: Contains images used for this tutorial

Description of system

Simulation trajectories contain a poly(ethylene terephthalate) (PET) trimer in acetone. PET is a common moiety in multilayer plastics that affords rigidity and would be useful to recycle as described in this link. Simulations of PET could capture different configurations that could be analyzed for solubility in different solvents using the COSMO-RS software package.

Visualizing the system (Optional)

To visualize the system, open VMD and run the command:

play "/Users/alex/tutorials/2020-1_Compute_energy_groups/vmd_plot_STRAP.sh"

(Note, your path may be different, marked in blue)

The image on the right shows the PET trimer in acetone; each monomer is colored in red, black, and blue. We would like to compute the LJ interactions for PET with itself and PET-acetone.

The image shows the output of "gmx energy" outputted in a typical Gromacs energy file (denoted as *.edr).

gmx energy -f "poly_solv_prod.edr" -s "poly_solv_prod.tpr" -o "energy.xvg"

By default, Gromacs does not output energies of specific groups. The procedure below will show you how to define energy groups, then compute energies between energy groups.

Step 1: Defining energy groups using gmx make_ndx

Run the following command in terminal to generate an index file containing PET trimer:

gmx make_ndx -f poly_solv_prod.tpr -o energy_group_index.ndx << INPUTS

keep 0

r PT3

r ACE

q

INPUTS

The command will create an index file with PET trimer (i.e. PT3) and acetone (ACE) as separate groups.


Step 2: Edit MDP file to include energy groups

The following function edits the production mdp file (prod.mdp) and adds specific energy groups; copy and paste this into terminal or include this as part of your ~/.bashrc

function mdp_add_energy_groups ()

{

input_mdp_file_="$1";

output_mdp_file_="$2";

energy_groups_="$3";

cp -r "${input_mdp_file_}" "${output_mdp_file_}";

echo "" >> "${output_mdp_file_}";

echo "; ENERGY GROUPS " >> "${output_mdp_file_}";

echo "energygrps = ${energy_groups_}" >> "${output_mdp_file_}";

echo "----- mdp_add_energy_groups -----";

echo "--> Created ${output_mdp_file_} from ${input_mdp_file_}";

echo "--> Added energy groups: ${energy_groups_}"

}

Then, add PET and ACE as an energy group by running:

mdp_add_energy_groups "prod.mdp" "prod_energy_groups.mdp" "PT3 ACE"

This function copies prod.mdp into a new file called prod_energy_groups.mdp and adds energy groups to the bottom:

; ENERGY GROUPS

energygrps = PT3 ACE

Step 3: Create a new *.tpr file with the new *.mdp file and rerun Gromacs

Gromacs needs a new *.tpr file with the energy groups, so create it using grompp:

gmx grompp -f "prod_energy_groups.mdp" -c "poly_solv_prod.gro" -p "poly_solv.top" -o "poly_solv_prod_energy_groups.tpr" -n "energy_group_index.ndx" -maxwarn 5

This should create a new file called "poly_solv_prod_energy_groups.tpr". Then, perform mdrun again using the "-rerun" flag:

gmx mdrun -v -nt 1 -s poly_solv_prod_energy_groups.tpr -rerun poly_solv_prod.xtc -e poly_solv_prod_energy_groups.edr

This should only take a couple of seconds to run, which will go through the trajectory and compute energies defined as "energygrps."

Step 4: Re-run "gmx energy" to get output of energy groups

gmx energy -f "poly_solv_prod_energy_groups.edr" -s "poly_solv_prod_energy_groups.tpr" -o "energy_groups.xvg"

Now, you have new energy groups:

  • "LJ-SR:PT3-PT3" - Lennard-Jones short-ranged interactions for PT3-PT3

  • "LJ-14:PT3-PT3" - Lennard-Jones 1,4-interactions for PT3-PT3

Ditto for ACE-ACE and PT3-ACE interactions. There are also coulomb terms (e.g. "Coul-SR"), but since the simulation was performed with Particle Mesh Ewald (PME) summation to quantify coulombic interactions, they cannot be decoupled as described in this Link.

Selecting "LJ-SR:PT3-PT3" will output LJ energies as a function of time, which might be useful to quantify things like solute-solvent interactions or quantify changes to LJ energies when modifying forcefields.


Now, you should be able to generate energy groups and re-compute energies with the "-rerun" command in Gromacs. This tool is useful for computing energies between solute-solvents and for checking changes in energies when modifying forcefield parameters.


Best of luck computing energies in Gromacs!

Last updated: 09/01/2021