Computing relative free energies with FEP using Q for multiple ligands

The script (prep_fepq.bash) exists with an accompanying manual that can be used to compute relative free energies with Q. This script is simple to use with a single ligand pair, but if one wants to do this for several ligands the setup can become rather tedious. Therefore, additional scripts were created that add another level of black-box on top of the prep_fepq.bash script. This guide describes the steps necessary to setup the calculations.

Samuel Genheden, June-July, 2011

Prerequisite

In a folder (I usually call it Prep) you must have the following

A pdb file for each protein-ligand complex, with hydrogen atoms added
An Amber prep file with the ending .in for each ligand

Before we start, make sure that the $QPREP path is set. Locally it should be,

export QPREP=/away/bio/Q/Q/Utils

Creating preparatory files

If the complex contains ions or co-factors, or if the ligand itself needs additional parameters, the following steps need to be executed.

First, we will create a file called leap_lines that tells leap to load the additional information. This will be added to the leap input file that the script is created. Here is an example for a Ca2+ ion.

cat << EOF > leap_lines loadamberprep ca.in loadamberparams ca.dat EOF
Next, we need to create a qlib file for the ion, which is the corresponding file in Q.

python $QPREP/in2qlib.py ca.in > ca.qlib
And then we need to create a qprep_lines file that tells Qprep5 to read this qlib-file

cat << EOF > qprep_lines readlib ca.qlib EOF

Creating a library of ligand pairs

The next step is to create a library of ligand pairs, one for each mutation that should be considered. It is important the first ligand in a pair have more or equal number of atoms than the second ligand, and that the atoms that should not be perturbed have the same name in both ligands.

First, we need to create a file called pairs. This file should contain one line for each pair. Each line should contain the residue names of the two ligands, separated by a dash (-). E.g.:

cat << EOF > pairs L01-L02 L01-L03 L02-L04 EOF

Next, execute

bash $QPREP/make_preplib.bash pairs

This script will take each pair and find out the atoms that differ between the two ligands. You need the specify the corresponding atoms, and if the atom should be perturbed to a dummy atom it can be left blank. The list of corresponding atoms for each pair will be stored in a file called prep_lib

Use pair library and execute prep_fepq.bash

Now we have all the preparatory files that are necessary to execute prep_fepq.bash. Therefore, type:

bash $QPREP/run_preplib.bash

The script will ask for the following information

The name of the protein
The desired size of the system
The equilibration length in ps
The production length in ps
The number of transformations (1 or 2)
The lambda-values at which the simulations should be done
The cluster where the calculations are to be executed

After these questions, the script will execute prep_fepq.bash for each of the ligand pairs. With the information from the questions and the prepartory steps above, it is possible to run the setup for each pair automatically. The script will create a folder called ../L01tl02, if the first ligand is called L01 and the second L02 and execute prep_fepq.bash there. When it is finnished will all the pair, the files in ../*/Fep_sim needs to be copied to the cluster.