top of page

How to Calculate the Free Energy of Methane in Water Using Gromacs with Cloud HPC

A step-by-step tutorial about using Gromacs to calculate the free energy on a cloud-based HPC platform with an example of methane in water.



Summary


This tutorial will guide the user to calculate a simple free energy change: the decoupling of van der Waals interactions between neutral methane and a box of water. The calculation uses the 2021-fosscuda-2019b version of GROMACS. In this tutorial, we will only focus on using GROMACS on Cloudam and calculate the specific practical steps of free energy, in order to reproduce the results of the free energy calculation for a very simple system.

Unlike the thermodynamic integration methods used in some other tutorials, this tutorial uses the gmx bar module of GROMACS to analyze the data, which was introduced in GROMACS version 4.5 (formerly known as g_bar). There are many practical applications for free energy calculations, some common ones include the solvation or hydration energy, and the binding free energy between a small molecule and some large acceptor biomolecule (usually a protein). Both of these processes involve add (introduction/coupling) or removal (decoupling/elimination) of the studied small molecule from the system and the calculation of the corresponding free energy change. In the free energy

calculations, two types of non-bonding interactions can be converted, Coulomb interactions and van der Waals interactions. GROMACS can also handle bond interactions, but for safety purposes, there will be not discussed.

In this tutorial, we will calculate the free energy of a quite simple process: turning off the Lennard-Jones interaction between the atomic sites of the molecule to be studied in water (refer to methane).


Step 1: Examine the topology file


1. Download the coordinate file of the system to be studied along with the topology file, and examine the topology file, where you will find the following points.


; Topology for methane in TIP3P
#include "oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
Methane 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB
chargeB massB
1 opls_138 1 ALAB CB 1 0.000 12.011
2 opls_140 1 ALAB HB1 2 0.000 1.008
3 opls_140 1 ALAB HB2 3 0.000 1.008
4 opls_140 1 ALAB HB3 4 0.000 1.008
5 opls_140 1 ALAB HB4 5 0.000 1.008
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
[ angles ]
; ai aj ak funct c0 c1 c2 c3
2 1 3 1
2 1 4 1
2 1 5 1
3 1 4 1
3 1 5 1
4 1 5 1
; water topology
#include "oplsaa.ff/tip3p.itp"
[ system ]
; Name
Methane in water
[ molecules ]
; Compound #mols
Methane 1
SOL 596

It can be noticed that all of the charges are set to zero. There is a reason for this setting: normally before the van der Waals interactions are turned off, the interaction of charges between solute and water is turned off first. If only the interaction of charges is retained and the Lennard-Jones term is turned off, the distance between the positive and negative charges approaching each other may be infinitely small, leading to a quite unstable system and even collapse possibly.


2. After examining the topology files, upload the two files to Cloudam for backup: create the Methane folder in General Workspace, and upload the two files into that folder, and finally rename the files to methane_water.gro and topol.top


Step 2: Develop a workflow


The workflow used is as follows:

1. Maximum speed down energy minimization

2. L-BFGS energy minimization

3. NVT balancing

4. NPT balancing

5. MD simulation

Two energy minimization steps are used here for better energy minimization of the initial structure. L-BFGS usually converges in advance and gives unstable systems, but satisfactory results are obtained when coupled with the most rapid descent method.


Step 3: Prepare mdp file


1. Create a new em_steep.mdp file for maximum speed down energy minimization with the following contents

Create a new em_l-bfgs.mdp file for L-BFGS energy minimization with the following contents:


; Run control
integrator = steep
nsteps = 5000
; EM criteria and other stuff
emtol = 100
emstep = 0.01
niter = 20
nbfgscorr = 10
; Output control
nstlog = 1
nstenergy = 1
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 1
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
; Temperature and pressure coupling are off during EM
tcoupl = no
pcoupl = no
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each
simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and
lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in
this case)
sc-power = 1
sc-sigma = 0.3
couple-moltype = Methane ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; No velocities during EM
gen_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12

2. Create a new em_l-bfgs.mdp file for L-BFGS energy minimization with the following contents


; Run control
integrator = l-bfgs
nsteps = 5000
define = -DFLEXIBLE
; EM criteria and other stuff
emtol = 100
emstep = 0.01
niter = 20
nbfgscorr = 10
; Output control
nstlog = 1
nstenergy = 1
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 1
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
; Temperature and pressure coupling are off during EM
tcoupl = no
pcoupl = no
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each
simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and
lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in
this case)
sc-power = 1
sc-sigma = 0.3
couple-moltype = Methane ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; No velocities during EM
gen_vel = no
; options for bonds
constraints = none ; L-BFGS doesn't work with constraints
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12

3. Create a new nvt.mdp file for NVT balancing with the following contents


; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 50000 ; 100 ps
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = system
tau_t = 1.0
ref_t = 300
; Pressure coupling is off for NVT
Pcoupl = No
tau_p = 0.5
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each
simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and
lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in
this case)
sc-power = 1
sc-sigma = 0.3
couple-moltype = Methane ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; Generate velocities to start
gen_vel = yes
gen_temp = 300
gen_seed = -1
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12

4. Create a new npt.mdp file for NPT balancing with the following contents


; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 50000 ; 100 ps
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = system
tau_t = 1.0
ref_t = 300
; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau_p = 1.0
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each
simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and
lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in
this case)
sc-power = 1
sc-sigma = 0.3
couple-moltype = Methane ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; Do not generate velocities
gen_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NVT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12

5. Create a new md.mdp file for MD simulation with the following contents


; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 500000 ; 1 ns
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = system
tau_t = 1.0
ref_t = 300
; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau_p = 1.0
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each
simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and
lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in
this case)
sc-power = 1
sc-sigma = 0.3
couple-moltype = Methane ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; Do not generate velocities
gen_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12

6. Upload these 5 files to Cloudam for backup: Create the MDP folder in General Workspace and upload these 5 files to that folder.


Step 4: Scripting the job


As mentioned before, the energy minimization will be performed in two phases: the first phase uses the fastest descent method and the second phase uses the L-BFGS method. This combination is able to get a fully energy-minimized structure, which is well-suited as a starting point for equilibration and for the subsequent data collection. Afterward, two stages will be equilibrated, the first one is an equal volume equilibrium (NVT) and the second one is an equal pressure equilibrium (NPT). Once the system is balanced, we can start the finished simulation.

1. The job.sh script reads as follows:


#!/bin/bash
module add GROMACS/2021.5-foss-2021b-CUDA-11.4.1-PLUMED-2.8.0
# Set some environment variables
FREE_ENERGY=/home/cloudam
echo "Free energy home directory set to $FREE_ENERGY"
MDP=$FREE_ENERGY/MDP
echo ".mdp files are stored in $MDP"
# A new directory will be created for each value of lambda and
# at each step in the workflow for maximum organization.
mkdir Lambda
cd Lambda
#################################
# ENERGY MINIMIZATION 1: STEEP #
#################################
echo "Starting minimization for lambda ..."
mkdir EM_1
cd EM_1
# Iterative calls to grompp and mdrun to run the simulations
gmx grompp -f $MDP/em_steep.mdp -c $FREE_ENERGY/Methane/methane_water.gro -p
$FREE_ENERGY/Methane/topol.top -o min.tpr
gmx mdrun -nt 2 -deffnm min
sleep 10
#################################
# ENERGY MINIMIZATION 2: L-BFGS #
#################################
cd ../
mkdir EM_2
cd EM_2
# We use -maxwarn 1 here because grompp incorrectly complains about use of a
plain cutoff; this is a minor issue
# that will be fixed in a future version of Gromacs
gmx grompp -f $MDP/em_l-bfgs.mdp -c ../EM_1/min.gro -p
$FREE_ENERGY/Methane/topol.top -o min.tpr -maxwarn 1
# Run L-BFGS in serial (cannot be run in parallel)
gmx mdrun -nt 1 -deffnm min
echo "Minimization complete."
sleep 10
#####################
# NVT EQUILIBRATION #
#####################
echo "Starting constant volume equilibration..."
cd ../
mkdir NVT
cd NVT
gmx grompp -f $MDP/nvt.mdp -c ../EM_2/min.gro -p $FREE_ENERGY/Methane/topol.top
-o nvt.tpr
gmx mdrun -nt 2 -deffnm nvt
echo "Constant volume equilibration complete."
sleep 10
#####################
# NPT EQUILIBRATION #
#####################
echo "Starting constant pressure equilibration..."
cd ../
mkdir NPT
cd NPT
gmx grompp -f $MDP/npt.mdp -c ../NVT/nvt.gro -p $FREE_ENERGY/Methane/topol.top -
t ../NVT/nvt.cpt -o npt.tpr
gmx mdrun -nt 2 -deffnm npt
echo "Constant pressure equilibration complete."
sleep 10
#################
# PRODUCTION MD #
#################
echo "Starting production MD simulation..."
cd ../
mkdir Production_MD
cd Production_MD
gmx grompp -f $MDP/md.mdp -c ../NPT/npt.gro -p $FREE_ENERGY/Methane/topol.top -t
../NPT/npt.cpt -o md.tpr
gmx mdrun -nt 2 -deffnm md
echo "Production MD complete."
# End
echo "Ending. Job completed for lambda"

2. Upload job.sh to Cloudam for backup: upload to jobs folder in General Workspace.


Step 5: Submit the GROMACS job


Start the management node in General Workspace and connect to it.



Use the cd jobs command to enter the /home/cloudam/jobs directory.



Use sbatch -N 1 -p g-t4-1 -c 4 job.sh to submit the GROMACS job (free energy calculation of methane in water), which is expected to take about 15 minutes.


Step 6: Analysis of results


In the working directory, simply collect all the md*.xvg files generated by mdrun, then run gmx bar :


gmx bar -f md*.xvg -o -oi -oh

The program prints out a lot of useful information on the screen and generates three output files: bar.xvg , barint.xvg and histogram.xvg. The screen output of gmx bar looks like the following:


Detailed results in kT (see help for explanation):
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 0.07 0.00 0.03 0.00 0.03 0.00 0.25 0.00
1 2 0.06 0.00 0.03 0.00 0.04 0.00 0.26 0.00
2 3 0.05 0.00 0.03 0.00 0.04 0.00 0.27 0.00
3 4 0.03 0.00 0.03 0.00 0.04 0.00 0.28 0.00
4 5 0.02 0.00 0.04 0.00 0.05 0.00 0.29 0.00
5 6 -0.00 0.01 0.04 0.00 0.05 0.00 0.30 0.01
6 7 -0.02 0.01 0.05 0.01 0.06 0.01 0.32 0.00
7 8 -0.06 0.00 0.06 0.00 0.07 0.00 0.35 0.00
8 9 -0.10 0.01 0.06 0.00 0.08 0.01 0.38 0.00
9 10 -0.15 0.01 0.08 0.01 0.10 0.01 0.41 0.00
10 11 -0.23 0.01 0.10 0.01 0.13 0.01 0.47 0.00
11 12 -0.32 0.01 0.10 0.01 0.14 0.01 0.51 0.01
12 13 -0.44 0.01 0.17 0.01 0.20 0.01 0.58 0.01
13 14 -0.57 0.01 0.15 0.01 0.17 0.01 0.57 0.00
14 15 -0.60 0.02 0.11 0.01 0.11 0.01 0.47 0.00
15 16 -0.53 0.01 0.06 0.01 0.06 0.01 0.35 0.00
16 17 -0.41 0.00 0.03 0.00 0.03 0.00 0.25 0.00
17 18 -0.28 0.00 0.02 0.00 0.02 0.00 0.18 0.00
18 19 -0.16 0.00 0.01 0.00 0.01 0.00 0.13 0.00
19 20 -0.05 0.00 0.00 0.00 0.00 0.00 0.10 0.00
Final results in kJ/mol:
point 0 - 1, DG 0.19 +/- 0.01
point 1 - 2, DG 0.16 +/- 0.00
point 2 - 3, DG 0.12 +/- 0.00
point 3 - 4, DG 0.09 +/- 0.00
point 4 - 5, DG 0.04 +/- 0.00
point 5 - 6, DG -0.00 +/- 0.01
point 6 - 7, DG -0.06 +/- 0.01
point 7 - 8, DG -0.15 +/- 0.01
point 8 - 9, DG -0.26 +/- 0.02
point 9 - 10, DG -0.39 +/- 0.02
point 10 - 11, DG -0.58 +/- 0.02
point 11 - 12, DG -0.79 +/- 0.02
point 12 - 13, DG -1.11 +/- 0.02
point 13 - 14, DG -1.43 +/- 0.02
point 14 - 15, DG -1.51 +/- 0.04
point 15 - 16, DG -1.33 +/- 0.02
point 16 - 17, DG -1.02 +/- 0.01
point 17 - 18, DG -0.70 +/- 0.01
point 18 - 19, DG -0.40 +/- 0.00
point 19 - 20, DG -0.13 +/- 0.00
total 0 - 20, DG -9.25 +/- 0.15

Conclusion


We have calculated a simple transition process: free energy change for decoupling of van der Waals interactions between simple solute molecules (uncharged methane) and the solvent (water), which has been calculated previously with high accuracy. Hope that this tutorial will give you a hand when dealing with complex systems.

Have fun with your simulations!


About Cloudam


Cloudam HPC is a one-stop HPC platform with 300+ pre-installed to deploy immediately. The system can smartly schedule compute nodes and dynamically schedule the software licenses, optimizing workflow and boosting efficiency for engineers and researchers in Life Sciences, AI/ML, CAE/CFD Simulations, Universities/Colleges, etc.


Partnered with AWS, Azure, Google Cloud, Oracle Cloud, etc., Cloudam powers your R&D with massive cloud resources without queuing.


You can submit jobs by intuitive templates, SLURM, and Windows/Linux workstations. Whether you are a beginner or a professional, you can always find it handy to run and manage your job.


There is a $30 Free Trial for every new user. Why not register and boost your R&D NOW?



bottom of page