Step 5: QM/MM → MM free energies

In step 3, you calculated the average difference in energy between the QM/MM and MM models using Monte Carlo sampling of the QM/MM model. In step 4, you saw how you can use λ to scale from the QM/MM model to the MM model. In this step, we will combine steps 3 and 4 together so that you can use Monte Carlo to sample the average energy difference at different values of λ. This average energy difference is useful, as it is actually equal to the gradient of the free energy with respect to λ. If you want to see why, take a look here.

If you calculate the gradients of the free energy across λ, and then integrate them, you will return the potential of mean force across λ, from which you can obtain the difference in free energy between the QM/MM and MM models of each side-chain. The program quantomm, is built using Sire and does exactly this. Quantomm performs Multiscale Monte Carlo simulations at eight different values of λ, calculating the average difference between the QM/MM and MM models for each value of λ. These are then integrated across λ to return the PMF, and thus the difference in free energy between the QM/MM and MM models. If you are interested, you can take a look at the quantomm python script.

In this step, you will use quantomm to calculate the difference in free energy between the MM and your chosen QM/MM model for both side-chains. Each quantomm simulation must be run in its own directory, so to start, create a directory for the PHE side-chain, and then change into that directory;

This will create a file called “freenrg.output” that will contain the free energies estimates from each of the five blocks of sampling, together with the potential of mean force (PMF) and final relative free energy. The final free energy is just the value of the PMF at λ = 1. The PMF is calculated together with an error, allowing the PMF and its “Maximum” and “Minimum” values to be calculated and printed. Your simulation is very short, and so the error is difficult to calculate. It is likely that the “Maximum” and “Minimum” values of the PMF are “nan”, indicating that the error could not be calculated.

EXERCISE 4 - Performing and analysing quantomm free energy simulations

Repeat the procedure above to perform a quantomm simulation for the TRP sidechain. Remember to create and change into a new directory for the simulation. You should use the same “config” file as for the PHE simulation above, and also use the same semi-empirical method as you used for PHE.

Analyse the output of both simulations using sire.app/bin/analyse_freenrg. Plot the PMFs for PHE and TRP on the same graph (ignoring the “Maximum” and “Minimum” values that are likely to be “nan”). Calculate the difference in free energy between the QM/MM and MM models for PHE and TRP (the values of the PMFs at λ = 1).

QUESTION 4.1 - How does the PMF for each base calculated in this step compare to the graphs of energy versus λ calculated in the last step (step 4)?

QUESTION 4.2 - How does the difference in free energy between your QM/MM and MM model calculated in this step compare to the average energy differences for that model calculated in step 3?

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

Now take a look at the help that comes with quantomm. To do this, type;

You can see from this that we need to supply the name of the residue that identifies the solute that will be modelled using QM/MM (in this case, “-l PHE”). You also need to supply the name of the topology and coordinate files that contain the solute in the water box (in this case, “-c ../PHE.crd” and “-t ../PHE.top”). You also need to supply the QM method to use to model the solute, e.g. “-q AM1” to use AM1, “-q PM6” to use PM6 etc. For our simulations, we will also supply our own config file that is called “config” and in the directory above (using “-C ../config”. Please take a look at this file to see what it does), and we will use the option “--intermolecular-only” to tell quantomm to use a QM model only for the intermolecular interaction between the side-chain and the water molecules. The complete command to run the simulation is therefore;

where you should replace “METHOD” with the semi-empirical method that you used in the last step. For example, I used PM6 in the last step, so I should type;

Running this simulation will take a long time (about 10 minutes) so you may want to take a break or check out some of the background information provided with this website.

When you run the simulation, you will see that quantomm prints out information about the moves that are being performed and the λ values that are being sampled. You should be able to see that quantomm is using eight λ values; 0.0, 0.142, 0.285, 0.429, 0.571, 0.714, 0.857 and 1.0. You should also see that quantomm is using four types of Monte Carlo move to sample each λ value;

(1)RigidBodyMC( maximumTranslation() = 0.1875 A, maximumRotation() = 2.14698 degrees ) - this performs random rigid body translation and rotation moves on the solute (PHE), with a maximum translation of 0.1875 angstroms, and maximum rotation of 2.14698 degrees.
(2)InternalMove() - this performs random bond, angle and dihedral moves on the solute, thereby sampling the intramolecular degrees of freedom of PHE. This uses the “flexibility” information that is automatically derived for PHE, and is printed further up the output (under the line “Flexibility of solute MolNum(1) equals”).
(3)RigidBodyMC( maximumTranslation() = 0.15 A, maximumRotation() = 15 degrees ) - this performs random rigid body translation and rotation moves on the solvent (water) molecules, with a maximum translation of 0.15 angstroms and maximum rotation of 15 degrees.
(4)VolumeMove( maximumVolumeChange() = 150.1 A^3 ) - this performs random changes in the size of the periodic box, with a maximum change in box volume of 150.1 cubic angstroms. This keeps the pressure constant, therefore we are sampling from the NPT (isothermal-isobaric) ensemble.

The quantomm simulation is run as a series of five blocks, with each block containing five multiscale Monte Carlo moves. The output of the quantomm simulation is written to the “output” directory. This contains a summary of the free energy results for each block (the files “output/results_XXXX.log”, where “XXXX” is the block number). This also contains PDB files for the coordinates of the atoms at λ = 0 and λ = 1 for each block (the files “output/coords_XXX_YYY.pdb”, where “XXX” is the block number and “YYY” is the λ value). The most important file is “output/freenrgs.s3”, which contains all of the average energy differences from every block and every value of λ. To combine these together to give the PMF across λ, type;