How to run a RESP calculation


Most information here is now obsolete: RESP fits are now done in an automatic way by antechanber. Follow this procedure.
  1. Optimise the molecule with QM methods, typically B3LYP/6-31G*.
    Avogadro can automatically add protons to a new molecule if they are missing.

  2. Calculate the electrostatic potential. For the Amber-99SB force field, this should be done at the HF/6-31G* level (B3LYP/6-31G* if you have metal sites). Other force fields require other levels of theory:
    HF/6-31G* for ff94, ff99, ff14SB, ff19SB, GAFF.
    B3L YP/cc-pVTZ for ff02
    B3L YP/cc-pVTZ with IEFPCM and a dielectric constant of 4 for ff03


    Example of a Gaussian input file:
    %Chk=btn.chk
    #P B3LYP/6-31g* Opt

    Biotin, HF, 6-31G*, 1/9-05

       -1    1
    c      33.34700000000000    7.55100000000000   14.58000000000000
    ...
    h      36.15300000000000   11.00700000000000   16.90000000000000

    --Link1--
    %Chk=btn.chk
    #P HF/6-31G*     SCF=Tight Geom=AllCheck Guess=Read
       Pop=MK IOp(6/33=2,6/41=10,6/42=17)

  3. Run antechamber on the output file of this calculation:

    antechamber -fi gout -fo prepi -c resp -i file.out -o res.in -rn RES -at amber -pf y

    where, you should insert the name of the output file from the previous point, and the residue name, (three letters) twice, first in lower-case and second in capital letters.
    Use -at amber, if want to use only standard Amber atom types. The default is to  get GAFF atom types, which is normally better
    -pf y indicates that intermediate files should be deleted.

  4. Often Amber gives you strange atom names. Then, you can use the program fixatomin to replace the wrong atom names. Insert the correct atom names into the file NEWPDB.PDB file and then read it into the program together with the NEWPDB.PDB file and the *.in file from antechamber.

  5. However, for the polarisable Amber ff02 force field, you still need to use this home page (because you need to remove the self-polarisation of the molecule). The procedure is described in Cieplak, Caldwell & Kollman, J. Comput. Chem. 22 (2001) 1048.

  6. If antechamber fails, you may also need to use the old version by hand.
    Alternatively, the program MCPB.py in AmberTools might be tested.

How to determine RESP charges for the polarisable Amber ff02 force field.

  1. Run a Gaussian geometry optimisation, either at the
    HF/6-31G* level or at the B3LYP/DZpdf/6-31G* level (with transition metals).

    Then run a single-point energy calculations with the B3LYP/cc-pVTZ level and
    include the following keywords in the command line.
    SCF=Tight Pop=MK IOp(6/33=2,6/41=10,6/42=17)

    If you already have a converged wavefunction at the appropriate level of theory, you can simply run:
    Guess=(Read,Only) Pop=(Minimal,MK) IOp(6/33=2,6/41=10,6/42=17)

    If the system includes transition metal, you also have to specify:
    Pop=ReadRadii
    and give a row for each metal at the end of the script, (with one blank line before and after) e.g.
    Cu 2.0

    IOp(6/33=2) makes Gaussian write out the potential points and potentials (do not change)
    IOp(6/41=10) specifies that 10 concentric layers of points are used for each atom (do not change)
    IOp(6/42) givest the density of points in each layer. A value of 17 gives about 2500 points/atom. Lower values may be needed for large molecules, since the programs cannot normally handle more than 100 000 potential points. A value of 10 gives about 1000 points/atom.

    Example of a full Gaussian script:
    %Chk=btn.chk
    #P HF/6-31g* Opt

    Biotin, HF, 6-31G*, 1/9-05

       -1    1
    c      33.34700000000000    7.55100000000000   14.58000000000000
    c      33.06400000000000    8.81200000000000   15.35700000000000
    ...
    h      36.15300000000000   11.00700000000000   16.90000000000000

    --Link1--
    %Chk=btn1.chk
    #P B3LYP/cc-pVTZ     SCF=Tight Geom=AllCheck Guess=Read
       Pop=MK IOp(6/33=2,6/41=10,6/42=17)

  2. Use the program changepot to convert the Gaussian output file to resp input files. All defaults are valid, except the name of the Gaussian file.
    This gives the files resp.in, resp.pot, and qin.

  3. Set qwt=0.0000 in the resp.in file.

    Then run resp:
    resp -O -i resp.in -o resp.out -q qin -e resp.pot

    This gives the first charges (in file qout). Check that they are identical to the input charges in qin (can be seen in resp.out).

  4. Run
    antechamber -i file.out -fi gout -o res.in -fo prepi -c resp -rn RES -at amber
    to get a *.in file for the molecule.
    You also get a PDB file, NEWPDB.PDB

  5. Run tleap to get the prmtop and prmcrd files for the molecule alone:
    tleap -s -f leaprc.ff02

    loadAmberPrep bt1.in
    loadamberparams btn.dat (optional)
    x=loadpdb NEWPDB.PDB
    check x
    saveamberparmpol x prmtop prmcrd
    quit
     
  6. Convert the prmtop file to the old format:
    new2oldparm <prmtop >p.for

  7. If amber reorders the atoms, you need to construct a file "order" with:
    #atoms
    1 line for each atom with number in the gaussian=resp files and in the amber=*.in file (arbitrary order)

    You can get a pdb file from the gaussian calculation by
    gausstopdb

    Otherwise, you do
    ln qout qout1

  8. Run the (undocumented program) changeQ_parm to move the charges from resp into the prmtop file:
    $AMBERHOME/src/resp/changeQ_parm

    It reads the file p.for and qout1 and writes the file p_mod.for, which is a prmtop file with charges taken from qout1.
    If the program is not compiled, you can compile it by make changeQ_parm in $AMBERHOME/src/resp.

  9. Run sander with the following sander.in file:
    Title
     &cntrl
      irest=0,ntx=1,
      imin=1,maxcyc=1,
      ntc=1,ntf=1,
      cut=999.0,
      ntpr=100,ntwx=0,ntwv=0,ntwe=0,
      ipol=1,iesp=1,
      igb=0,ntb=0,
      scnb=2.0,scee=1.2
     &end
     &ewald
      indmeth=1
     &end

    The (undocumented) parameter iesp=1 forces the program to read the file esp.dat and write out the two files esp.induced and esp.qm-induced. The latter is the ESP(QM)-ESP(induced dipole), which we need for the new Resp calculations. Thus run:

    ln resp.pot esp.dat
    sander -O -i sander.in -o sander.out -c prmcrd -p p_mod.for

  10. Set up some links to make the iterations working
    cp qout1 q0

    1. Now you can iterate in the following way (in each subsequent iteration, you should increase the highlighted number after q with 1):

      resp -O -i resp.in -o resp.out -q qout -e esp.qm-induced
      reorder_qout
      cp qout1 q1
      less resp.out

      \mv p_mod.for p.for
      $AMBERHOME/src/resp/changeQ_parm
      sander -O -i sander.in -o sander.out -c prmcrd -p p_mod.for

      You should iterate (repete these lines) until the charges (qin and qout) in resp.out are identical. Alternatively, you may check it by
      diff q7 q8

      You can do it automatically by this script:

      cp qout1 q0
      xold=0
      for x in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
      do
               resp -O -i resp.in -o resp.out -q qout -e esp.qm-induced
         reorder_qout
         cp qout1 q"$x"
         diff q"$xold" q"$x" >script.tmp
         if [ ! -s script.tmp ] ; then exit ; fi
         cat script.tmp
         \      mv p_mod.for p.for
         $AMBERHOME/src/resp/changeQ_parm
         sander -O -i sander.in -o sander.out -c prmcrd -p p_mod.for
         xold=$x
      done

  11. Control if there are any polar hydrogens bound to the same atoms or if there are identical polar heavy atoms (e.g. in carboxylic groups).
    If so, constrain these to have the same charge (in resp.in1).

    Example (    water)

    Water
     &cntrl  ihfree=0, qwt=0.0005, iqopt=2
     &end
     1.0
    Mol1
        0    3
        8    0
        1    0
        1    2

    The lines contain the following information:
    Title
    Control keywords
    Weight of the molecule
    Name of the molecule
    Total charge, number of atoms
    Atom number  Fitting_number

    Valid fitting_numbers are:
    -n  fixed charge (to the one in file qin)
     0  free charge
     n  the charge is constrained to be the same as the one of atom #n.

    The keywords mean:
    ihfree (0 = hydrogen atoms are restrained; 1 = not restrained)
    qwt (wheight of the constraint. Use 0.0005 for stage 1 and 0.0001 for stage 2)
    iqopt (1 = set all initial charges to 0; 2 = read initial charges from -q (qin)
    irstrnt (0 = harmonic restraints; 1 = hyperbolic restraints, default)
    nmol (number of molecules, i.e. different conformations)
     
  12. Run resp:
    resp -O -i resp.in1 -o resp.out1 -q qout -e esp.qm-induced
     
  13. Fix the charge of all atoms except carbons with more than one hydrogen (methyl and methylene groups).
         Constrain hydrogens bound to the same carbon to have equal charges.
         Change qwt=0.001.
    Do this in resp.in2.
    Example:

    Methanol
     &cntrl  ihfree=0, qwt=0.001, iqopt=2
     &end
     1.0
    Mol1
        0    6
        6    0
        1    0
        1    2
        1    2
        8   -1
        1   -1
  14. Run resp again:
    resp -O -i resp.in2 -o resp.out2 -q qout -e esp.qm-induced
  15. The results are found in file resp.out2
    Insert these charges into the *.in file created by antechanber above.
    Note that antechanber renumbers the atoms, so you need to check that you use the right atom, but the file order gives the right order.

    You can check the charges by running
    changeparm
    p_mod.for
    p
    pdb
    m
    prmcrd
    q

How to determine RESP charges - old version by hand (normally obsolete)

1. Optimise the geometry, either at the


2. For Turbomole, you set $esp_fit resp in the control file and then run
ridft-resp -proper >logm

For Gaussian:
Run a Gaussian single-point energy calculation, either at the
HF/6-31G* level (without transition metals) or at the
B3LYP/DZpdf/6-31G* level.

Include the following keywords in the command line.
SCF=Tight Pop=MK IOp(6/33=2,6/41=10,6/42=17)

If you already have a converged wavefunction at the appropriate level of theory, you can simply run:
Guess=(Read,Only) Pop=(Minimal,MK) IOp(6/33=2,6/41=10,6/42=17)

If the system includes transition metal, you also have to specify:
Pop=ReadRadii
and give a row for each metal at the end of the script, (with one blank line before and after) e.g.
Cu 2.0

IOp(6/33=2) makes Gaussian write out the potential points and potentials (do not change)
IOp(6/41=10) specifies that 10 concentric layers of points are used for each atom (do not change)
IOp(6/42) givest the density of points in each layer. A value of 17 gives about 2500 points/atom. Lower values may be needed for large molecules, since the programs cannot handle more than 100 000 potential points. A value of 10 gives about 1000 points/atom.

Example of a full Gaussian script:
%Chk=g.chk
#P B3LYP/Gen Geom=AllCheck Guess=(Read,Only)
   Pop=(Minimal,MK,ReadRadii)
   IOp(6/33=2,6/41=10,6/42=17)

! Gen basis set
@/teo/ulf/Basis/stdbas/N

! Non-default atomic radii
 Zn 2.0
 Cu 2.0
 

3. Use the program changepot to convert the output file to resp input files. All defaults are valid, except the name of the Gaussian file.
This gives the files resp.in, resp.in1, resp.pot, and qin.
 

4. Control if there are any polar hydrogens bound to the same atoms or any symmetry related other atoms.
If so, constrain these to have the same charge (in resp.in), e.g. for water:

Water
 &cntrl  ihfree=0, qwt=0.0005, iqopt=2
 &end
 1.0
Mol1
    0    3
    8    0
    1    0
    1    2

The lines contain the following information:
Title
Control keywords
Weight of the molecule
Name of the molecule
Total charge, number of atoms
Atom number  Fitting_number

Valid fitting_numbers are:
-n  fixed charge (to the one in file qin)
 0  free charge
 n  the charge is constrained to be the same as the one of atom #n.

The keywords mean:
ihfree (0 = hydrogen atoms are restrained; 1 = not restrained)
qwt (wheight of the constraint. Use 0.0005 for stage 1 and 0.0001 for stage 2)
iqopt (1 = set all initial charges to 0; 2 = read initial charges from -q (qin)
irstrnt (0 = harmonic restraints; 1 = hyperbolic restraints, default)
nmol (number of molecules, i.e. different conformations)
 

5. Run resp:
resp -O -i resp.in -o resp.out -q qin -e resp.pot
 

6. Fix the charge of all atoms except carbons with more than one hydrogen (methyl and methylene groups).
     Constrain hydrogens bound to the same carbon to have equal charges.
     Change qwt=0.001.
Do this in resp.in1.
There is now already a template resp.in1 file (from changepot) that is correct, provided that hydrogen atoms come after their respective heavy atom. Othewise, you have to this by hand.
cp qout qin1

Example:

Methanol
 &cntrl  ihfree=0, qwt=0.001, iqopt=2
 &end
 1.0
Mol1
    0    6
    6    0
    1    0
    1    2
    1    2
    8   -1
    1   -1

7. Run resp again:
resp -O -i resp.in1 -o resp.out1 -q qin1 -e resp.pot

The result is found in file resp.out1

Add extra charge constraints (e.g. of capping ACE and NME groups):
Add these directly after the atoms line, without any blank line between.
Example with four charge constraints for a Zn(His)2(H2O)2 site, with capped His groups (constraint in bold at the end):

Charges for HHWW site from 1VLI, UR 23/5-11                                    
 &cntrl  ihfree=0, qwt=0.0005, iqopt=2
 &end
 1.0
Mol1
    2   65
    1    0
    6    0
    1    0
    1    0
    6    0
    8    0
    7    0
    1    0
    6    0
    1    0
    6    0
    1    0
    1    0
    6    0
    7    0
    1    0
    6    0
    1    0
    7    0
    6    0
    1    0
    6    0
    8    0
    7    0
    1    0
    6    0
    1    0
    1    0
    1    0
    1    0
    6    0
    1    0
    1    0
    6    0
    8    0
    7    0
    1    0
    6    0
    1    0
    6    0
    1    0
    1    0
    6    0
    7    0
    1    0
    6    0
    1    0
    7    0
    6    0
    1    0
    6    0
    8    0
    7    0
    1    0
    6    0
    1    0
    1    0
    1    0
   30    0
    8    0
    1    0
    1   61
    8    0
    1    0
    1   64
    6   0.0
    1    1    1    2    1    3    1    4    1    5    1    6
    6   0.0
    1   24    1   25    1   26    1   27    1   28    1   29
    6   0.0
    1   30    1   31    1   32    1   33    1   34    1   35
    6   0.0
    1   53    1   54    1   55    1   56    1   57    1   58




How to run a Resp calculation with multiple conformations
  1. Run changepot on all Gaussian files for each configuration.

    2. Save the potentials in different files
     
  2. Concatenate all potential files to one single file

    3. cat files >espot   (if you use this file name, you do not need to have a -e flag)
     
  3. Add

    4. nmol=number_of_configurations
    to the  &cntrl list in resp.in
     
  4. Add a section:

    5.  1.0
    Mol1
        0    3
        8    0
        1    0
        1    2

    for each molecule.
    The first line is a general weight for the molecule (preferably the Boltzmann weight based on the total energy (in solution or in a MD simulation).
    The second line is the name of the molecule.
    The third row is the total charge and the number of atoms as usual.
    The rest of the rows are the normal atom number - charge constraints rows.
    There should be one blank row between each molecule.
     

  5. Finally, after two blank rows, add a section which equates the charge of the same atom in all molecules. 
    The program /teo/ulf/Bin/respconfig generates these row (fixed format).
    You just give the number of atoms and number of molecules.
    The format is two rows for each atom (16i5), the first with the number of constraints = number of molecules, and the second is the constraint: molecule_1 atom_1, molecule_2 atom_2, ...
        6

    6.     1    1    2    1    3    1    4    1    5    1    6    1
        6
        1    2    2    2    3    2    4    2    5    2    6    2
        6
        1    3    2    3    3    3    4    3    5    3    6    3
    ...

  6. Then run the two stage RESP as usual.
    It seems that there are some error in the code if the weights differ from 1.0. Then the residual sum can become negative, giving NAN for the stddev and rrms. I do not know if the results are correct.

RESP with Turbomole

  1. Optimise the structure with QM.

  2. If it is an UHF calculation, you need to have natural orbitals. If you do not:
    Insert into the control file the following keywords and rerun ridft (dscf):
    $natural orbitals  file=natural
    $natural orbital occupation  file=natural

  3. Run makepoints to get the points where you will calculate the electrostatic potential (ESP).
    Select Turbomole (t), Ball (b), 3.0, 0.3, ChelpG (c), 10000, 8 Å, and Turbomole (non-default values highlighted).
    This gives the file esppoints.
    Note that moloch only writes out at the most 10 000 ESP points. If you want more than that, you need to construct several files and concatenate the results.
    Therefore, you should probably modify the distance between the points (in green above), until you get approximately 10 000 points.

    Old version (<5.8):

  4. Insert the suggested rows into the control file (in the beginning to avoid already existing keywords).

  5. Run moloch (RHF) or moloch2 (UHF) to get the potentials and also the Mulliken charges.
    moloch2 is stupid and does not copy the file esppoints to the new directory. Therefore, you have to do:
    moloch2
    cd moloch_n_dir
    cp ../esppoints .
    moloch2
    The calculation takes ~30 minutes and the result is in moloch_n_dir/moloch.out.

  6. This file can be read by changepot as detailed above.
With the new (>5.7) version of Turbomole, some changes are needed:
  1. Insert in the control file
    $pointval geo=point pot fmt=xyz
     
    and then the list of coordinates from the file esppoints
  2. Run
    ridft -proper >logm
  3. The results is in the file tp.xyz

How to ensure that RESP can read more than 9999 potential points
by Matteo Guglielmi 090824
vim ${AMBERHOME}/src/antechamber/espgen.c
# 207c207
# <       fprintf(fpout, "%5d%5d%5d\n", tmpint1, tmpint2, 0);
# ---
# >       fprintf(fpout, "%10d%10d%10d\n", tmpint1, tmpint2, 0);

vim ${AMBERHOME}/src/resp/resp.f
# 927,928c927,928
# <         read(10,'(2i5)') inat,nesp
# <         WRITE(6,'(/,t2,a,i3,/,t2,a,i5,/,t2,a,i5)')
# ---
# >         read(10,'(2i10)') inat,nesp
# >         WRITE(6,'(/,t2,a,i3,/,t2,a,i5,/,t2,a,i10)')

vim ${AMBERHOME}/src/sander/sander.f
# 1046,1048c1046,1048
# <    read (dat_unit,'(3i5)')inat,nesp,idum
# <    write(6,'(t2,''inat = '',i5)')inat
# <    write(6,'(t2,''nesp = '',i5)')nesp
# ---
# >    read (dat_unit,'(3i10)')inat,nesp,idum
# >    write(6,'(t2,''inat = '',i10)')inat
# >    write(6,'(t2,''nesp = '',i10)')nesp
# 1050,1051c1050,1051
# <    write(new_unit,'(2i5)')inat,nesp
# <    write(minus_new_unit,'(2i5)')inat,nesp
# ---
# >    write(new_unit,'(2i10)')inat,nesp
# >    write(minus_new_unit,'(2i10)')inat,nesp

vim ~/bin/CHANGEPOT/changepot.f
# 216c216
# <   111  format(2i5)
# ---
# >   111  format(2i10)


... and then recompile the whole amber9 package.

References
1. CI Bayly et al, J. Phys. Chem. 97(93)10269.
2. WD Cornell et al, J. Am. Chem. Soc. 115(93)9620.