Generating parameters using RESP charges from GAMESS-US

This tutorial is for anyone wanting to generate restricted electrostatic potential (RESP) charges for a custom ligand or other molecule. The starting point is an AMBER .prmtop and .inpcrd file for the molecule of interest. Steps on how to generate these can be found in this excellent tutorial by Ross Walker.

You can download a tar file containing all the input and expected output for this tutorial here. To extract them, just run:

wget http://www-wales.ch.cam.ac.uk/chris/resp_tutorial.tar.bz
tar xvfj resp_tutorial.tar.bz

The process we're going to use is very simple:

1. Run AMBGMIN with initial charges (GROUPROTATION moves)
2. Calculate RESP charges for global minimum using GAMESS-US
3. Run AMBGMIN with new set of RESP charges
4. Re-calculate RESP charges for new global minimum
5. Repeat as many times as you like!

Before you start

Here is a quick list of things to check now - without all of these, you will not be able to proceed! If you don't understand any of these steps, please ask someone for help now rather than later!

• you must have the latest version of the group SVN repository SCRIPTS directory:

cd ~/svn/SCRIPTS
svn update

• you must be working on a machine with GAMESS-US compiled! It should be in ~csw34/gamess, you can check with ls:

csw34@volkhan:~> ls ~csw34/gamess/rungms
/home/csw34/gamess/rungms

At the moment, GAMESS-US is compiled on mek-quake, clust and volkhan - NOT sinister.

• you must have an empty ~/scr directory. If it doesn't exist, make it!

mkdir ~/scr

• you must have your $AMBERHOME variable set - i.e. in ~/.bashrc:

export AMBERHOME=/home/csw34/amber11

• you must have $AMBERHOME/exe in your $PATH - i.e. in ~/.bashrc

export PATH=$PATH:$AMBERHOME/exe

• you must have AMBGMIN compiled and in your $PATH (you should be able to run it from anywhere by typing 'AMBGMIN')

Setting up

To set up, you need a few files and to make a directory:

1. make build_structure directory:

mkdir build_structure

This directory should contain your .prmtop and .inpcrd files for the molecule you want charges for.

2. AMBGMIN input files

You need a few files so that GMIN can run happily:

• data - the GMIN data file. Feel free to adjust as you wish, just make sure it still samples a lot of conformations during the run (make sure the Markov energy is changing).

COMMENT DEBUG
SLOPPYCONV 1.0D-4
TIGHTCONV 1.0D-5
MAXERISE 1.0D-4
DUMPINT 100
UPDATES 1000
MAXIT 20000 60000
MAXBFGS 1.0
CENTRE
TRACKDATA
GROUPROTATION 1
DUMPSTRUCTURES
TEMPERATURE 1.5
SAVE 100
STEPS 5000 1.0
STEP 0.0
RADIUS 300
AMBER9 coords.inpcrd inpcrd

• min.in - the AMBER input file specifying the potential. You might want to change the value of igb here, it's up to you. Take a look at the Sander input section in the AMBER manual here for info on each of the options. Note that ifswitch is not a standard AMBER keyword and so you won't find it in the manual!

STOP
 &cntrl
  imin   = 1,
  ncyc = 1,
  maxcyc = 1,
  igb = 2, saltcon=0.1,
  ntb    = 0,
  cut    = 999.99,
  ifswitch = 1,
 /

• atomgroups - the group definition file for the GROUPROTATION moves in GMIN. You'll need to make this yourself, as described under GROUPROTATION in the documentation here.

GROUP COO 9 10 2 1.0 0.6
11
12
GROUP NH2 14 15 2 1.0 0.6
16
17
GROUP AMIDE 18 19 7 1.0 0.4
20
21
22
23
26
27
28
GROUP CH2CH31 3 2 6 1.0 0.4
1
32
33
34
35
36
GROUP CH2CH32 3 4 6 1.0 0.4
5
38
39
40
41
42
GROUP 2CH2CH3 6 3 15 1.0 0.3
1
2
4
5
32
33
34
35
36
37
38
39
40
41
42
GROUP O2CH2CH3 7 6 16 1.0 0.3
1
2
3
4
5
32
33
34
35
36
37
38
39
40
41
42

• gmin_job - qsub file to start AMBGMIN jobs

#PBS -q s1

cd $PBS_O_WORKDIR
AMBGMIN > amber.out

Ok! With all those files, you should be good to go! For the tutorial, simply extract the zip as explained at the top - all these files should be placed correctly for you.

Calculating the charges

From here on, I'm going to assume you're using the files in resp_tutorial.tar.bz linked from the top. This means that the initial input files are OTV.prmtop and OTV.inpcrd. If your files are called something different, just replace OTV with whatever you have in all of the below! Don't forget to change the total charge too!

It's very easy to start the iterative calculation:

~/svn/SCRIPTS/AMBER/iterative_resp/control_script.sh OTV -1 3

This will start a run for:

• a system called OTV (build_structure contains OTV.prmtop and OTV.inpcrd)
• which has a total charge of -1 (from COO- in this case)
• for 3 iterations

You can take a look at the GMIN runs in iteration_X/gmin, and the GAMESS-US jobs in iteration_X/gamess. It's worth taking a look at the end of the output files (GMIN_out and OTV.log) in these directories to make sure they're running ok!

This could take a while to run, go get some coffee/tea/water :)

Important output files

Once the script has finished running, the important output files are:

• OTV_final.prep - the AMBER prep file containing the final RESP charges. You can use this file to create AMBER input for complexes containing the new molecule using tleap. Ask for help if you're not sure how to do this!
• iteration_3/gmin/coords.prmtop - an AMBER topology file containing the final RESP charges. This is useful if you'd like to run the molecule by itself. Note that if you run more than 3 iterations, you should look in the final directory.