Up to this point, you have been calculating energy differences between the QM/MM and MM models of the side-chains. What we want, by the end of this practical, is to calculate the free energy difference between the QM/MM and MM models. This will allow us to ‘correct’ the relative hydration free energy calculated using the MM models to give us the relative hydration free energy of the side-chains using the QM/MM models.

There are many ways of calculating free energy differences. The method we will use in this practical is called thermodynamic integration (described in detail here). To calculate the free energy difference between the QM/MM and MM models of the side-chains, we need to construct a reaction coordinate that maps between these two representations. This reaction coordinate is called λ, and is defined via the equation;

E(λ) = ( 1 - λ ) EQM/MM + λ EMM

At λ = 0, E(λ) = EQM/MM. At λ = 1, E(λ) = EMM. By varying λ between 0 and 1, we morph between the QM/MM and MM models of the side-chains.

This is illustrated in the simple python script, lambda.py. If you are curious, you can take a look at this script. This script uses the above equation to calculate the energy of one of the side-chains using a specified semi-empirical model, at a specified value of λ. You use the script by typing;

where “SIDECHAIN” is the name of the side-chain, “METHOD” is the name of the semi-empirical method you want to use, and “LAMBDA_VALUE” is the value of λ (between 0 and 1). For example, to calculate the energy of the tryptophan side-chain using PM6 at λ = 0.5, type;

$SIRE/bin/python lambda.py SIDECHAIN METHOD LAMBDA_VALUE

EXERCISE 3 - Performing an energy scan across λ

Using the same semi-empirical method you chose for the last exercise, calculate the value of E(λ) for each of the side-chains at the following values of λ; 0.0, 0.2, 0.4, 0.6, 0.8, 1.0. Place these values into a spreadsheet. Calculate the energy relative to that for λ = 0 (to do this, subtract E(λ=0) from each E(λ) value). Plot the relative energy for each side-chain versus λ on the same graph.

If you get stuck or need help, you can take a peek at this model answer for exercise 3.

$SIRE/bin/python lambda.py TRP PM6 0.5

← Previous                        Next →