Tutorial

Using the SFED python script with AMBER

This tutorial is assuming you're using Ambertools 18, you might need to modify some scripts to run in later version if anything has changed.

1) Choose a PDB file

Choose a PDB file for any structure. In this case we’ll use ethene. OpenBabel can be used to convert between different file formats if you can’t find a PDB file.

Create a directory

$ mkdir ethene

Then create a second directory in ethene in which you’ll parametrise the molecule

$ mkdir 1-prm

And place your PDB file in there

$ cp ethene.pdb ethene/1-prm/

2) Parametrise the PDB file

Create a file called runleap.in to write LEap script – this will parametrise the pdb.

runleap.in :

source leaprc.gaff
MOL = loadmol2 ethene.mol2
check MOL
loadamberparams  ethene.frcmod
saveoff MOL ethene.lib
saveamberparm MOL ethene.prmtop  ethene.rst7
quit

Then we’ll write a bash script to run antechamber and parmchk2 before using tleap in order to generate a mol2 file and assign charges, atom types and check if the relevant forcefield parameters are present.

prmtise.sh:

antechamber -i ethene.pdb -fi pdb -o ethene.mol2 -fo mol2 -c bcc -s 2
parmchk2 -i  ethene.mol2 -f mol2 -o ethene.frcmod
tleap -f runleap.in

Run the prmtise.sh script in the terminal using:

$ sh prmtise.sh

3) 1D-RISM - Solvent Preparation

Now the solvent distrbution around ethene must be computed. 1D-RISM is used to compute the solvent susceptibility which is then used as an input in 3D-RISM.

Make a directory in the ethene\ folder:

$ mkdir 2-1D-RISM

You'll be preparing the solvent here.

We'll use the SPC water model with Na+ and Cl- ions at concentrations 55.5M, 0.005M and 0.005M respectively at 298K using the Kovalenko-Hirate (KH) closure.

DRISM is used to ensure the dielectric constant is enforced throughout the calculation.

SPC_NaCl.inp:

&PARAMETERS  
OUTLST='x', THEORY='DRISM', CLOSUR='KH',  
!grid  
NR=16384, DR=0.025, routup=384, toutup=0,  
!MDIIS  
NIS=20, DELVV=0.3, TOLVV=1.e-12,  
!iter  
KSAVE=-1, KSHOW=1, maxste=10000,  
!ElStat  
SMEAR=1, ADBCOR=0.5,  
!bulk solvent properties  
TEMPER=298, DIEps=78.497,  
NSP=3  
/  
&SPECIES
DENSITY=55.5d0,  
MODEL="$AMBERHOME/dat/rism1d/model/SPC.mdl"  
/  
&SPECIES  
DENSITY=0.005d0,  
MODEL="$AMBERHOME/dat/rism1d/model/Na+.mdl"  
/  
&SPECIES  
DENSITY=0.005d0,  
MODEL="$AMBERHOME/dat/rism1d/model/Cl-.mdl"  
/

Then run the following command in the terminal: $ rism1d SPC_NaCl (not $ rism1d SPC_NaCl.inp - notice the lack of file extension) You should then see a file called SPC_NaCl.xvv - this is the solvent susceptibility.

4) 3D-RISM

Finally, we run the 3D-RISM package. We need to point it to the geometry, parameter and solvent files. We then need to output total, direct and pair correlation functions as well as the potential as .dx files.

Firstly, navigate back to the ethene/ folder and make a new directory:

$ mkdir 4-3D-RISM

And then write a 3D-RISM script (if relevant, this is the command you want to set up to queue and run on your cluster):

runrism3d.sh:

rism3d.snglpnt --pdb ../1-prm/ethene.pdb --prmtop ../1-prm/ethene.prmtop --closure KH --guv ethene_G --huv ethene_H --cuv ethene_C --uuv ethene_U --gf --pc+ --xvv ../2-1D-RISM/SPC_NaCl.xvv > ethene.out

Notice we tag the outputs with G, H, C, U. This becomes relevant when using the SFED script.

5) SFED

Navigate back to the ethene\ folder and make a new folder:

$ mkdir 5-SFED

Make sure you have the SFED code downloaded and placed here. Then run:

$ python sfed.py ../4-3D-RISM -i ethene -o ethene_SFED -c KH -p 0.03342285869 -T 298 You should then see ethene_SFED.dx in your folder - this is the solvation free energy density in OpenDX format.