MM/PBSA

Using mmpbsa.py

To run MMPBSA.Py you need:
prmtop files for complex, receptor, ligand
coordinate file from the MD simulation
input file for the MMGBSA calculations

<< mmgbsa.in >>
 &general
   verbose = 2 
   keep_files = 1 
   endframe = 9999
   interval = 1 
   startframe = 1 
/
 &gb
   surften = 0.00542
   surfoff = 0.92
   saltcon = 0.0 
   igb = 2

In the queue script this line is for the MMPBSA calculation, depending on the cluster change accordingly.

srun MMPBSA.py.MPI -i mmgbsa.in -y coordinate_file(in the name use * for a wildcard) -cp complex.prmtop  -rp receptor.prmtop -lp ligigand.prmtop -eo output.csv

The good thing about this script is that it finds automatically atoms belonging to complex, receptor and ligand. So there is no need to mess with any other files. And this script can do a lot more and it is straightforward to set it in the input file.  

The manual for this can be found at http://ambermd.org/tutorials/advanced/tutorial3/py_script/files/MMPBSA_Python_Manual.pdf

Paulius 23/10-15

Instructions how to run the MM/PBSA method with Amber 8.0 to obtain a ligand-binding affinity.
The method is described in Kollman, Acc. Chem. Res. 33 (2000) 889, Kuhn & Kollman J. Med. Chem. 43 (2000) 3786, and Kuhn, J. Med. Chem. 48 (2005) 4040.
  1. Start with a crystal structure with a ligand bound to the protein.
    Possibly, you need to mutate the ligand (by hand or with Spartan).
  2. Obtain parameters for the ligand using Antechamber.
    For ff-03, you should first optimise the ligand with the HF/6-31G** method and then calculate the ESP using B3LYP-IEFPCM/cc-pVTZ method (sample G98 input file).
    antechamber -i file.out -fi gout -o res.in -fo prepi -c resp -rn RES -at amber
    Here, you should insert the file name, and the residue name, (three letters) twice, first in lower-case and second in capital letters. -at amber, indicates that you will use only standard Amber atom types (otherwise you will get gaff atom types).
  3. Run leap for the ligand+protein to obtain prmtop and prmcrd files for the protein solvated in an octahedral box (tleap sample input file).
  4. Run 4 MD simulations (see sample input files below):
    1. 100 steps of steepest decent MM minimisation with all heavy atoms restrained to
      the original QM/MM structure.
    2. 20 ps MD simulation at constant pressure, still with heavy atoms restrained.
    3. 50 ps MD equilibration with constant pressure
    4. 200 ps MD equilibration with constant volume
    5. 200 ps data collection with constant volume
      sander -i sander.in1 -o sander.out1 -r mdrest1 -c prmcrd -ref prmcrd
      sander -i sander.in2 -o sander.out2 -r mdrest2 -c mdrest1 -ref prmcrd
      sander -i sander.in3 -o sander.out3 -r mdrest3 -c mdrest2
      sander -i sander.in4 -o sander.out4 -r mdrest4 -c mdrest3
      sander -i sander.in5 -o sander.out5 -r mdrest5 -c mdrest4 -e mden5 -x mdcrd5
      sander sample input files

  5. Ensure that the ligand is at the centre of the periodic box in all simulations (replace 252 below with the residue number of the ligand):
    mv mdcrd5 mdcrd5_unwrap
    cat >ptraj.in <<EOF
    trajin mdcrd5_unwrap
    trajout mdcrd5
    image origin center familiar com :252
    go
    EOF
    ptraj prmtop ptraj.in
    Check that ptraj recognizes the mask (Amber9 often does not; use Amber8 then).

  6. Set up the files needed for the first MM/PBSA run (base these on a pdb file from the last snapshot:
    changeparm <<EOL
    prmtop
    p
    pdb5
    c
    mdcrd5
    20
    q
    EOL
    mkdir E
    cd E
    grep -v WAT ../pdb5 >pdb5
    and check that leap does not add any new atoms):
    1. prmtop file for complex without water (modify the leap file and strip water from pdb file), compl.prmtop
    2. prmtop file for protein without water (modify the leap file and strip water and ligand from pdb file), prot.prmtop
    3. prmtop file for ligand alone (modify the leap file and grep the ligand from pdb file), lig.prmtop

  7. Construct a charge and a radius file, with the default names delphi.crg and delphi.siz.
    You can construct these files by running changeparm on the compl.prmtop file, command del.
    Select the residue numbe of the ligand, parse radii, some appropriate coordinates (e.g. compl.prmcrd).
    [Already done automatic:
    In these two files, you should finally delete the residue number of the ligand (so that the file works also for the ligand alone, when it has residue number 1), i.e. change
      O3    BT1 495    1.400
    to
      O3    BT1        1.400]
  8. Run the first MM/PBSA calculation (everything except the normal-mode calculation) using these structures. You have to modifie the input file slightly (mm_pbsa.in1 sample file).
    mm_pbsa.pl mm_pbsa.in1 >mm_pbsa.out1
    This may gives rise to various errors. See below for some solutions.

  9. Calculate the entopies:

    NB This  method for calculating the entropy term to ligand binding leads in many cases to a rather large standard deviation in the results for the binding affinity. The standard devitation may be reduced by replacing this point 8 with the method described below.

    Set up the files needed for the normal-mode calculation.
    You need to truncate the protein to contain only residues within 8 Å of the ligand.
    This is done with changepdb command fix on the file with the complex without water molecules.
    Ensure that you do it residue-wise and write out the coordinates.
    Unfortunately, this probably gives problem at the termini of the protein (leap insists to put charged terminals to all proteins). The only way to avoid this seems to be to edit the leaprc file and remove the corresponding addPdbResMap commands. Moreover, you have to remove bonds between non consecutive residues. These can be found by running leap, command check x and remove the long bonds by deleteBond. Base these always on protonated pdb files from the original calculation (from leap or changeparm) and check that leap does not add any new atoms.
    1. changepdb
      compl.pdb
      fix
      temp
      m
      n
      y


      name of the ligand
      residue number
      8
      n
      y
      compl_tr.pdb
      y
      q
    2. Construct prot_tr.pdb by removing the ligand.
    3. cp /home/bio/AMBER/Amber8/dat/leap/cmd/leaprc.ff03 .
    4. Remove all the rows under addPdbResMap, or only the few rows actually used.
    5. Insert TER into the compl_tr.pdb file
    6. tleap -s
      source leaprc.ff03
      loadAmberPrep bt1.in
      loadamberparams btn.dat
      x=loadpdb compl_tr.pdb
      deleteBond x.16.16 x.17.1
      check x
      saveamberparm x compl_tr.prmtop prmcrd
      quit
    7. Repeate for protein.
    8. Run MM/PBSA on these files with the properly modified mm_pdbsa.in2 file (note that the atom selection in the coordinate construction is probably much more complicated now, but it can be copied from the temp file)
      mm_pbsa.pl mm_pbsa.in2 >mm_pbsa.out2
    9. Analyse all the results.

New method to calculate the entropy

Recently, we published a new method to increase the precision of the entropy term by a factor of ~3 (so that it no longer limits the precision of MM/PBSA) by using a fixed buffer region. Thereby, the questionable use of a distant-dependent dielectric constant is also avoided.

The method is described in:
J. Kongsted & U. Ryde (2009) "An improved method to predict the entropy term with the MM/PBSA approach", J. Comput. Aided Mol Design, 23, 63-71; DOI: 10.1007/s10822-008-9238-z.

Set up the files needed for the normal-mode calculation.
  1. First, construct a pdb file (pdb5) with the protein complex AND water molecules:
    changeparm
    prmtop
    p
    pdb5
    c
    mdcrd5
    20
    q

  2. Then, truncate the protein to contain only residues within 12 Å of the ligand. Ensure that you do it residue-wise and write out the coordinates.
    changepdb
    pdb5
    fix
    temp
    m
    n
    y


    name of the ligand
    residue number, if necessary
    12
    n
    y
    compl_tr.pdb
    y
    q


  3. Construct prot_tr.pdb by removing the ligand.
    grep -v BT compl_tr.pdb > prot_tr.pdb

  4. cp /home/bio/AMBER/Amber9/dat/leap/cmd/leaprc.ff03 .

  5. Remove all the lines in leaprc.ff03 under addPdbResMap, or only the few rows actually used.
    This is because you otherwise get charged terminal groups where the protein is cut  (leap insists to put charged terminals to all proteins).

  6. Insert TER into the pdb file:
    changepdb <<EOF
    compl.pdb
    ter
    1
    w

    q
    EOF

    No, not needed:
    Remove bonds between non consecutive residues. These can be found by running leap, command check x and remove the long bonds by deleteBond.
    If you put the warnings of the type
    WARNING: There is a bond of 4.807912 angstroms between:
    -------  .R<ILE 54>.A<C 18> and .R<ARG 55>.A<N 1>
    into a file, the program fixdelbon makes the deleteBond entries for you.

    deleteBond x.16.16 x.17.1

    tleap -s
    source leaprc.ff03
    loadAmberPrep bt1.in
    loadamberparams btn.dat
    x=loadpdb compl_tr.pdb
    check x
    saveamberparm x compl_tr.prmtop prmcrd
    quit

    Also, check that leap does not add any new atoms at loadpdb.

    Repeat after adding deleteBond.

  7. Repeat for prot_tr.pdb.
    Typically, you should delete the last deleteBond entry for the protein, compared to the complex file.
    tleap -s
    source leaprc.ff03
    loadAmberPrep bt1.in
    loadamberparams btn.dat
    x=loadpdb prot_tr.pdb
    deleteBond x.16.16 x.17.1
    saveamberparm x prot_tr.prmtop prmcrd
    quit

  8. Ensure that you already have the lig.prmtop file (from mm_pbsa.in1)

  9. Modify the mm_pdbsa.in2 file:
    Change NUMBER_REC_GROUPS, RSTART/RSTOP list from to the temp file, and put
    DIELEC = 1 since we now have a buffer around the active region.
    Finally, set GC = 0, because we will not use MAKECRD to build the coordinates

  10. Now, we should divide the system into an active and a buffer region. The active region (which is minimised) consists of the atoms within 8 Å from the ligand and the buffer (which is not minimised)  is simply the rest. However, all water molecules are considered as beloning to the buffer,even if they are within 8 Å of the ligand. Use changepdb  to truncate compl_tr.pdb within 8 Å from the ligand. Save the truncation information using create_input format in the file temp2.
    changepdb
    compl_tr.pdb
    ren

    fix
    temp2
    c
    n
    y


    BTN1
    8
    n
    y
    temp2.pdb
    y
    q

  11. cp /home/bio/Data/Amber/mm_pbsa_loc_createinput.pm  .
    and then add the belly from temp2 into this file, remember that water molecules should NOT be included in the list (already fixed by changepdb).
    The belly should be changed in four places (look for ATOM), the two first will all belly lines, and the two following without the ligand in the belly. Sample file.                                                                                                                     

  12. Run changecrd, command l, on the mdcrd file(s) to generate the *crd files.
    This is necessary, because the water molecules will change in each snapshot. Therefore, this program adds the coordinates of the right number of closest water molecule to the crd file of each snapshot.
    To run the program, you need to know the number of atoms in the MD simulations, and the number of the first water atom in those simulations:
    You have to run a separate calculation for each mdcrd file (if you have several), using different offsets.

    changecrd
    40370
    mdcrd5
    l
    mm_pbsa.in2
    7833
    0
    q

  13. Run MM/PBSA (note the slightly different program name):
    mm_pbsa_loc.pl mm_pbsa.in2 >mm_pbsa.out2

Changes needed in the code for this procedure


cp mm_pbsa.pl  mm_pbsa_loc.pl
(Note that in Amber8, mm_pbsa.pl in src and exe differ; use the one in exe)
Change:
26c26
< use mm_pbsa_loc_createinput  qw(create_input);
---
> use mm_pbsa_createinput  qw(create_input);

In nmode.f (around line 296):
         write(6,*) 'Thermo analysis not supported for belly calc.'

!    Code added by Jacob Kongsted 2007 for new MM/PBSA entropy
! JK
         nvecs =  3*natsys - 6
         call thermo(natsys,nvecs,ilevel,x(mx),x(mamass),x(mcval), &
               x(mh),x(mh+ns3),x(mh+2*ns3), &
               x(mh+3*ns3),t,patm)
!    End of new code

In make_crd_hg.f:
        parameter (NMO=400)

Sample of the mm_pbsa_createinput_example.pm file

Only start of the file.
Changes marked in yellow.

#
# Module with functions to create input for mm_pbsa "helper" programs
#
# Holger Gohlke
#   17.04.2002
#
########################################################################

package mm_pbsa_loc_createinput;
require Exporter;

@ISA         = qw(Exporter);
@EXPORT      = qw(create_input);
@EXPORT_OK   = qw();
$VERSION     = 1.00;

########################################################################

use strict;
use lib ("$ENV{AMBERHOME}/src/mm_pbsa");
use mm_pbsa_global qw(@GC_PAR);

########################################################################

sub create_input(){
###################

  # Parameters: \%GEN,\%DEC,\%SAN,\%DEL,\%GBO,\%MOL,\%NMO,\%MAK,\%PRO,\$TRA
  my $r_gen = shift;
  my $r_dec = shift;
  my $r_san = shift;
  my $r_del = shift;
  my $r_gbo = shift;
  my $r_mol = shift;
  my $r_nmo = shift;
  my $r_mak = shift;
  my $r_pro = shift;
  my $r_tra = shift;

  print "\n";
  print "=>> Creating input\n";

  if($r_gen->{"MM"} || $r_gen->{"GB"}){
    &create_sander_input($r_gen,$r_dec,$r_san,$r_del,$r_gbo);
  }

  if($r_gen->{"PB"}){
    if($r_del->{"PROC"} == 1){
      &create_delphi_input($r_del);
    }
    elsif($r_del->{"PROC"} == 2){
      &create_pbsa_input($r_gen, $r_san, $r_del);
    }
  }

  if($r_gen->{"GC"} || $r_gen->{"AS"}){
    &create_makecrd_input($r_gen,$r_mak,$r_tra);
  }

  if($r_gen->{"NM"}){
    &create_nmode_input($r_gen,$r_nmo);
  }
 
  return;
}

sub create_nmode_input(){
#########################

  # Parameters: \%GEN,\%NMO
  my $r_gen = shift;
  my $r_nmo = shift;

  print "    Nmode input\n";

  my ($dielc, $maxcyc, $drms, $nmdrms);
  $dielc = $r_nmo->{"DIELC"};
  $maxcyc = $r_nmo->{"MAXCYC"};
  $drms = $r_nmo->{"DRMS"};
  $nmdrms = $drms * 100.0; # Increased drms for nmode due to differences
                          # in gradient calc. with respect to sander.

  my $suf_list = "";
  $suf_list .= "_com " if($r_gen->{"COMPLEX"});
  my $suf;
  foreach $suf (split/ +/,$suf_list){
    my $name = "sanmin" . $suf . ".in";
    open(OUT,">$name") || die("    $name not opened\n");
    print OUT "File generated by mm_pbsa.pl\n";
    print OUT " &cntrl\n";
    print OUT "  ntxo   = 0,\n";
    print OUT "  ntf    = 1,       ntb    = 0, ibelly=1\n";
    print OUT "  dielc  = $dielc,\n";
    print OUT "  cut    = 99.0,    nsnb   = 99999,\n";
    print OUT "  scnb   = 2.0,     scee   = 1.2,\n";
    print OUT "  imin   = 1,       maxcyc = 2000,\n";
    print OUT "  ncyc   = 0,       drms   = $drms\n";
    print OUT " &end\n";
    print OUT " --Belly\n";
    print OUT " ATOM    25   132\n";
    print OUT " ATOM   180   373\n";
    print OUT " ATOM   389   399\n";
    print OUT " ATOM   460   479\n";
    print OUT " ATOM   494   690\n";
    print OUT " ATOM   753   891\n";
    print OUT " ATOM   954  1072\n";
    print OUT " ATOM  1216  1234\n";
    print OUT " ATOM  1247  1292\n";
    print OUT " ATOM  1386  1399\n";
    print OUT " ATOM  1471  1484\n";
    print OUT " ATOM  1582  1612\n";
    print OUT "END\n";
    print OUT "END\n";

    close(OUT);

    $name = "nmode" . $suf . ".in";
    open(OUT,">$name") || die("    $name not opened\n");
    print OUT "File generated by mm_pbsa.pl\n";
    print OUT " &data\n";
    print OUT "  ntx    = 0,\n";
    print OUT "  ntrun  = 1, ibelly=1,       nvect  = 0,\n";
    print OUT "  drms   = $nmdrms,\n";
    print OUT "  dielc  = $dielc,     idiel  = 1,\n";
    print OUT "  scnb   = 2.0,     scee   = 1.2\n";
    print OUT " &end\n";
    print OUT " --Belly\n";
    print OUT " ATOM    25   132\n";
    print OUT " ATOM   180   373\n";
    print OUT " ATOM   389   399\n";
    print OUT " ATOM   460   479\n";
    print OUT " ATOM   494   690\n";
    print OUT " ATOM   753   891\n";
    print OUT " ATOM   954  1072\n";
    print OUT " ATOM  1216  1234\n";
    print OUT " ATOM  1247  1292\n";
    print OUT " ATOM  1386  1399\n";
    print OUT " ATOM  1471  1484\n";
    print OUT " ATOM  1582  1612\n";
    print OUT "END\n";
    print OUT "END\n";
    close(OUT);
  }

# JK

  {
  my $suf_list = "";
  $suf_list .= "_rec " if($r_gen->{"RECEPTOR"});
  my $suf;
  foreach $suf (split/ +/,$suf_list){
    my $name = "sanmin" . $suf . ".in";
    open(OUT,">$name") || die("    $name not opened\n");
    print OUT "File generated by mm_pbsa.pl\n";
    print OUT " &cntrl\n";
    print OUT "  ntxo   = 0,\n";
    print OUT "  ntf    = 1,       ntb    = 0, ibelly=1\n";
    print OUT "  dielc  = $dielc,\n";
    print OUT "  cut    = 99.0,    nsnb   = 99999,\n";
    print OUT "  scnb   = 2.0,     scee   = 1.2,\n";
    print OUT "  imin   = 1,       maxcyc = 2000,\n";
    print OUT "  ncyc   = 0,       drms   = $drms\n";
    print OUT " &end\n";
    print OUT " --Belly\n";
    print OUT " ATOM    25   132\n";
    print OUT " ATOM   180   373\n";
    print OUT " ATOM   389   399\n";
    print OUT " ATOM   460   479\n";
    print OUT " ATOM   494   690\n";
    print OUT " ATOM   753   891\n";
    print OUT " ATOM   954  1072\n";
    print OUT " ATOM  1216  1234\n";
    print OUT " ATOM  1247  1292\n";
    print OUT " ATOM  1386  1399\n";
    print OUT " ATOM  1471  1484\n";
    print OUT "END\n";
    print OUT "END\n";

    close(OUT);

    $name = "nmode" . $suf . ".in";
    open(OUT,">$name") || die("    $name not opened\n");
    print OUT "File generated by mm_pbsa.pl\n";
    print OUT " &data\n";
    print OUT "  ntx    = 0,\n";
    print OUT "  ntrun  = 1, ibelly=1,       nvect  = 0,\n";
    print OUT "  drms   = $nmdrms,\n";
    print OUT "  dielc  = $dielc,     idiel  = 1,\n";
    print OUT "  scnb   = 2.0,     scee   = 1.2\n";
    print OUT " &end\n";
    print OUT " --Belly\n";
    print OUT " ATOM    25   132\n";
    print OUT " ATOM   180   373\n";
    print OUT " ATOM   389   399\n";
    print OUT " ATOM   460   479\n";
    print OUT " ATOM   494   690\n";
    print OUT " ATOM   753   891\n";
    print OUT " ATOM   954  1072\n";
    print OUT " ATOM  1216  1234\n";
    print OUT " ATOM  1247  1292\n";
    print OUT " ATOM  1386  1399\n";
    print OUT " ATOM  1471  1484\n";
    print OUT "END\n";
    print OUT "END\n";
    close(OUT);
  }
  }

Energy decomposition

With GB/PBSA, the interaction energy can be decomposed, both residue-wise and residue-pair-wise
a sample input file is shown below.
Note that only GB is allowed here, and even MAKECRD must be run in a separate calculation. Special parameters for GBSA must also be used.
Finally, it should be noted that all data in the DECOMP input must be ranges; thus LIGRES 1-1 must be given rather than 1.

The output file (*_statistics.out) is a very large file with four sections:
=>> COMPLEX
=>> RECEPTOR
=>> LIGAND
=>> DELTA

In each section, one line for each residue is given, and in each of these lines, 50 numbers are given: A consecutive number from 0, residue number, and then 24 energies and standard deviations in 2*3*8 sets, namely the eight normal values:
INT    - internal energy (should be 0)
VDW    - van der Waals energy
ELE    - electrostatic energy
GAS    - sum of these three values
GB     - the GB solvation energy
SUR    - the surface area energy (non-polar solvation energy)
GBSOL  - sum GB + SUR
GBTOT  - sum GAS+GBSOL

The three different components of each value: S, B, and T, are (for amino acids) contributions from the side-chain and back-bone atoms, as well as their total value.
In fact, there are twice the number of energies, but every second energy is zero and can be ignored.

    Number    Residue       SINT                  BINT                  TINT                  SVDW                  BVDW                  TVDW   SELE                  BELE                  TELE                  SGAS                  BGAS                  TGAS                   SGB  BGB                   TGB                SGBSUR                BGBSUR              TGBSUR                SGBSOL                BGBSOL                TGBSOL                SGBTOT                BGBTOT                TGBTOT




    G98 input file for charge calculation with ff03

    Note the empty line at the end of the input.

    %Chk=btn2.chk
    #P HF/6-31g** Opt

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

       -1    1
    c      33.34700000000000    7.55100000000000   14.58000000000000
    c      33.06400000000000    8.81200000000000   15.35700000000000
    c      34.33700000000000    9.27000000000000   16.14000000000000
    c      34.10000000000000   10.39700000000000   17.15500000000000
    c      35.40800000000000   11.01800000000000   17.70300000000000
    c      36.01200000000000   10.28500000000000   18.93300000000000
    c      36.51100000000000    9.73300000000000   21.32800000000000
    c      37.33800000000000   10.93700000000000   19.27100000000000
    c      37.59500000000000   10.64500000000000   20.76600000000000
    c      37.00500000000000   12.89700000000000   20.51000000000000
    o      34.63100000000000    7.35800000000000   14.25700000000000
    o      32.46300000000000    6.78000000000000   14.30900000000000
    n      36.73800000000000   14.11100000000000   20.81400000000000
    n      37.43200000000000   11.96400000000000   21.34000000000000
    n      36.91700000000000   12.34400000000000   19.28000000000000
    s      35.04900000000000   10.35900000000000   20.46900000000000
    h      32.75800000000000    9.59300000000000   14.65200000000000
    h      32.24600000000000    8.61200000000000   16.05800000000000
    h      35.07800000000000    9.61800000000000   15.41100000000000
    h      33.52900000000000    9.99100000000000   17.99800000000000
    h      35.19800000000000   12.05400000000000   17.99300000000000
    h      36.17600000000000    9.29400000000000   18.49500000000000
    h      36.76600000000000    8.69500000000000   21.08500000000000
    h      36.48100000000000    9.85600000000000   22.41700000000000
    h      38.17600000000000   10.67600000000000   18.61500000000000
    h      38.56600000000000   10.19400000000000   20.99800000000000
    h      36.59000000000000   12.89300000000000   18.49500000000000
    h      37.61300000000000   12.22800000000000   22.30000000000000
    h      34.73500000000000    8.40300000000000   16.68000000000000
    h      33.51800000000000   11.18700000000000   16.66700000000000
    h      36.15300000000000   11.00700000000000   16.90000000000000
    h      36.85000000000000   14.43000000000000   21.76800000000000

    --Link1--
    %Chk=btn2.chk
    #P B3LYP/cc-pvtz SCRF=(IEFPCM,Read) Geom=AllCheck Guess=Read
       Pop=(Minimal,MK) IOp(6/33=2,6/41=10,6/42=17)

    eps=4

    Leap sample input file

    source leaprc.ff03
    loadAmberPrep btn.in
    loadamberparams btn.dat
    x=loadpdb input.pdb
    bond x.4.SG x.83.SG
    solvateOct x TIP3PBOX 10
    saveamberparm x prmtop prmcrd
    quit

    Sander sample input files

    The sander.in1 file
    100 steps of steepest decent MM minimisation with all heavy atoms restrained to
    the original QM/MM structure.
    ---
    Title
     &cntrl
      irest=0,ntx=1,
      imin=1,maxcyc=100,drms=0.0001,ntmin=2,
      ntc=2,ntf=2,
      cut=8.0,
      ntpr=100,ntwx=0,ntwv=0,ntwe=0,
      ipol=0,igb=0,
      scnb=2.0,scee=1.2,
      ntr=1,restraint_wt=100,
      restraintmask="!:WA= & !@H="
     &end


    The sander.in2 file
    20 ps MD simulation at constant pressure, still with heavy atoms restrained.
    ---
    Title
     &cntrl
      irest=0,ntx=1,
      nstlim=10000,dt=0.002,
      temp0=300.0,ntt=1,tautp=1.0,
      ntc=2,ntf=2,
      cut=8.0,
      ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
      ntb=2,ntp=1,pres0=1.0,taup=1.0,
      ipol=0,igb=0,
      scnb=2.0,scee=1.2,
      ntr=1,restraint_wt=50,
      restraintmask="!:WA= & !@H="
     &end

    The sander.in3 file
    50 ps MD equilibration with constant pressure
    ---
    Title
     &cntrl
      irest=1,ntx=5,
      nstlim=25000,dt=0.002,
      temp0=300.0,ntt=1,tautp=1.0,
      ntc=2,ntf=2,
      cut=8.0,
      ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
      ntb=2,ntp=1,pres0=1.0,taup=1.0,
      ipol=0,igb=0,
      scnb=2.0,scee=1.2,
      ntr=0
     &end

    The sander.in4 file
    200 ps MD equilibration with constant volume
    ---
    Title
     &cntrl
      irest=1,ntx=5,
      nstlim=100000,dt=0.002,
      temp0=300.0,ntt=1,tautp=1.0,
      ntc=2,ntf=2,
      cut=8.0,
      ntpr=5000,ntwx=0,ntwv=0,ntwe=0,
      ntb=1,
      ipol=0,igb=0,
      scnb=2.0,scee=1.2,
      ntr=0
     &end

    The sander.in5 file
    200 ps data collection with constant volume
    ---
    Title
     &cntrl
      irest=1,ntx=5,
      nstlim=100000,dt=0.002,
      temp0=300.0,ntt=1,tautp=1.0,
      ntc=2,ntf=2,
      cut=8.0,
      ntpr=5000,ntwx=5000,ntwv=0,ntwe=5000,iwrap=1,
      ntb=1,
      ipol=0,igb=0,
      scnb=2.0,scee=1.2,
      ntr=0
     &end

    The first mm_pbsa.in1 sample file; Delphi version - Use this!

    It is essentially taken from /home/bio/AMBER/Amber8/src/mm_pbsa/Examples (but files merged and set up to have all files in the same directory).
    You should change issues marked in yellow.

    #
    # Input parameters for mm_pbsa.pl
    # This example just generates snapshots from a trajectory file
    #
    # Holger Gohlke
    # 08.01.2002
    #
    ################################################################################
    @GENERAL
    #
    # General parameters
    #   0: means NO; >0: means YES
    #
    #   mm_pbsa allows to calculate (absolute) free energies for one molecular
    #     species or a free energy difference according to:
    #
    #     Receptor + Ligand = Complex,
    #     DeltaG = G(Complex) - G(Receptor) - G(Ligand).
    #
    #   PREFIX - To the prefix, "{_com, _rec, _lig}.crd.Number" is added during
    #            generation of snapshots as well as during mm_pbsa calculations.
    #   PATH - Specifies the location where to store or get snapshots.
    #
    #   COMPLEX - Set to 1 if free energy difference is calculated.
    #   RECEPTOR - Set to 1 if either (absolute) free energy or free energy
    #              difference are calculated.
    #   LIGAND - Set to 1 if free energy difference is calculated.
    #
    #   COMPT - parmtop file for the complex (not necessary for option GC).
    #   RECPT - parmtop file for the receptor (not necessary for option GC).
    #   LIGPT - parmtop file for the ligand (not necessary for option GC).
    #
    #   GC - Snapshots are generated from trajectories (see below).
    #   AS - Residues are mutated during generation of snapshots from trajectories.
    #   DC - Decompose the free energies into individual contributions
    #        (only works with MM and GB).
    #
    #   MM - Calculation of gas phase energies using sander.
    #   GB - Calculation of desolvation free energies using the GB models in sander
    #        (see below).
    #   PB - Calculation of desolvation free energies using delphi (see below).
    #   MS - Calculation of nonpolar contributions to desolvation using molsurf
    #        (see below).
    #        If MS == 0, nonpolar contributions are calculated with the LCPO method
    #        in sander.
    #   NM - Calculation of entropies with nmode.
    #
    PREFIX                xxx
    PATH                  ./
    #
    COMPLEX               1
    RECEPTOR              1
    LIGAND                1
    #
    COMPT                 compl.prmtop
    RECPT                 prot.prmtop
    LIGPT                 lig.prmtop
    #
    GC                    1
    AS                    0
    DC                    0
    #
    MM                    1
    GB                    1
    PB                    1
    MS                    1
    #
    NM                    0
    #
    ################################################################################
    @MAKECRD
    #
    # The following parameters are passed to make_crd_hg, which extracts snapshots
    #   from trajectory files. (This section is only relevant if GC = 1 OR AS = 1 above.)
    #
    #   BOX - "YES" means that periodic boundary conditions were used during MD
    #         simulation and that box information has been printed in the
    #         trajecotry files; "NO" means opposite.
    #   NTOTAL - Total number of atoms per snapshot printed in the trajectory file
    #            (including water, ions, ...).
    #   NSTART - Start structure extraction from NSTART snapshot.
    #   NSTOP - Stop structure extraction at NSTOP snapshot.
    #   NFREQ - Every NFREQ structure will be extracted from the trajectory.
    #
    #   NUMBER_LIG_GROUPS - Number of subsequent LSTART/LSTOP combinations to
    #                       extract atoms belonging to the ligand.
    #   LSTART - Number of first ligand atom in the trajectory entry.
    #   LSTOP - Number of last ligand atom in the trajectory entry.
    #   NUMBER_REC_GROUPS - Number of subsequent RSTART/RSTOP combinations to
    #                       extract atoms belonging to the receptor.
    #   RSTART - Number of first receptor atom in the trajectory entry.
    #   RSTOP - Number of last receptor atom in the trajectory entry.
    #   Note: If only one molecular species is extracted, use only the receptor
    #         parameters (NUMBER_REC_GROUPS, RSTART, RSTOP).
    #
    BOX                   YES
    NTOTAL                36606
    NSTART                21
    NSTOP                 40
    NFREQ                 1
    #
    NUMBER_LIG_GROUPS     1
    LSTART                2622
    LSTOP                 3862
    NUMBER_REC_GROUPS     1
    RSTART                1
    RSTOP                 2621
    #
    #################################################################################
    @TRAJECTORY
    #
    # Trajectory names
    #
    #   The following trajectories are used to extract snapshots with "make_crd_hg":
    #   Each trajectory name must be preceeded by the TRAJECTORY card.
    #   Subsequent trajectories are considered together; trajectories may be
    #     in ascii as well as in .gz format.
    #   To be able to identify the title line, it must be identical in all files.
    #
    TRAJECTORY            ./mdcrd5
    #
    ################################################################################
    @PB
    #
    # PB parameters (this section is only relevant if PB = 1 above)
    #
    #   The following parameters are passed to the PB solver.
    #   Additional parameters (e.g. SALT) may be added here.
    #   For further details see the delphi and pbsa documentation.
    #
    #   PROC - Determines which method is used for solving the PB equation:
    #          If PROC = 1, the delphi program is applied. If PROC = 2,
    #           the pbsa program of the AMBER suite is used.
    #   REFE - Determines which reference state is taken for PB calc:
    #          If REFE = 0, reaction field energy is calculated with EXDI/INDI.
    #            Here, INDI must agree with DIELC from MM part.
    #          If REFE > 0 && INDI > 1.0, the difference of total energies for
    #            combinations EXDI,INDI and 1.0,INDI is calculated.
    #            The electrostatic contribution is NOT taken from sander here.
    #   INDI - Dielectric constant for the molecule.
    #   EXDI - Dielectric constant for the surrounding solvent.
    #   SCALE - Lattice spacing in no. of grids per Angstrom.
    #   LINIT - No. of iterations with linear PB equation.
    #   PRBRAD - Solvent probe radius in A (e.g. use 1.4 with the PARSE parameter set
    #     and 1.6 with the radii optimized by R. Luo)
    #
    #   Parameters for pbsa only
    #
    #   RADIOPT - Option to set up atomic avity radii for molecular surface calculation
    #     and dielectric assignment. A value of 0 uses the cavity radii from the prmtop file.
    #     A value of 1 sets up optimized cavity radii at the pbsa initialization phase.
    #     The latter radii are optimized for model compounds of proteins only; use cautions
    #     when applying these radii to nucleic acids.
    #
    #   Parameters for delphi only
    #
    #   FOCUS - If FOCUS > 0, subsequent (multiple) PERFIL and SCALE parameters are
    #     used for multiple delphi calculations using the focussing technique.
    #     The # of _focussing_ delphi calculations thus equals the value of FOCUS.
    #   PERFIL - Percentage of the lattice that the largest linear dimension of the
    #            molecule will fill.
    #   CHARGE - Name of the charge file.
    #   SIZE - Name of the size (radii) file.
    #
    #   SURFTEN / SURFOFF - Values used to compute the nonpolar contribution Gnp to
    #                  the desolvation according to Gnp = SURFTEN * SASA + SURFOFF.
    #
    #
    PROC                  1
    REFE                  0
    INDI                  1.0
    EXDI                  80.0
    LINIT                 500
    PRBRAD                1.4
    #
    FOCUS                 0
    PERFIL                90.0
    SCALE                 2.0
    CHARGE                delphi.crg
    SIZE                  delphi.siz

    #
    SURFTEN               0.00542    ! 0.0072 recommended for Amber charges; 0.005 for Parse charges recommended in Ras-Raf MM/PBSA tutorial
    SURFOFF               0.92
    #
    ################################################################################
    @MM
    #
    # MM parameters (this section is only relevant if MM = 1 above)
    #
    #   The following parameters are passed to sander.
    #   For further details see the sander documentation.
    #
    #   DIELC - Dielectricity constant for electrostatic interactions.
    #           Note: This is not related to GB calculations.
    #
    DIELC                 1.0
    #
    ################################################################################
    @GB
    #
    # GB parameters (this section is only relevant if GB = 1 above)
    #
    #   The first group of the following parameters are passed to sander.
    #   For further details see the sander documentation.
    #
    #   IGB - Switches between Tsui's GB (1), Onufriev's GB (2, 5).
    #   GBSA - Switches between LCPO (1) and ICOSA (2) method for SASA calc.
    #          Decomposition only works with ICOSA.
    #   SALTCON - Concentration (in M) of 1-1 mobile counterions in solution.
    #   EXTDIEL - Dielectricity constant for the solvent.
    #   INTDIEL - Dielectricity constant for the solute
    #
    #   SURFTEN / SURFOFF - Values used to compute the nonpolar contribution Gnp to
    #                   the desolvation according to Gnp = SURFTEN * SASA + SURFOFF.
    #
    IGB                   2
    GBSA                  1
    SALTCON               0.00
    EXTDIEL               80.0
    INTDIEL               1.0
    #
    SURFTEN               0.00542
    SURFOFF               0.92
    #
    ################################################################################
    @MS
    #
    # Molsurf parameters (this section is only relevant if MS = 1 above)
    #
    #   PROBE - Radius of the probe sphere used to calculate the SAS.
    #           Since Bondi radii are already augmented by 1.4A, PROBE should be 0.0
    #
    PROBE                 0.0
    #
    ################################################################################
    @PROGRAMS
    #
    # Program executables
    #
    # DELPHI                /sw/bio/Bin/delphi   on toto7; next line on docenten
    DELPHI                /sw/pkg/bio/Bin/delphi
    #
    ################################################################################

    The second mm_pbsa.in2 sample file

    For the normal-mode calculations.
    You should change issues marked in yellow.

    #
    # Input parameters for mm_pbsa.pl
    # This example just generates snapshots from a trajectory file
    #
    # Holger Gohlke
    # 08.01.2002
    #
    ################################################################################
    @GENERAL
    #
    # General parameters
    #   0: means NO; >0: means YES
    #
    #   mm_pbsa allows to calculate (absolute) free energies for one molecular
    #     species or a free energy difference according to:
    #
    #     Receptor + Ligand = Complex,
    #     DeltaG = G(Complex) - G(Receptor) - G(Ligand).
    #
    #   PREFIX - To the prefix, "{_com, _rec, _lig}.crd.Number" is added during
    #            generation of snapshots as well as during mm_pbsa calculations.
    #   PATH - Specifies the location where to store or get snapshots.
    #
    #   COMPLEX - Set to 1 if free energy difference is calculated.
    #   RECEPTOR - Set to 1 if either (absolute) free energy or free energy
    #              difference are calculated.
    #   LIGAND - Set to 1 if free energy difference is calculated.
    #
    #   COMPT - parmtop file for the complex (not necessary for option GC).
    #   RECPT - parmtop file for the receptor (not necessary for option GC).
    #   LIGPT - parmtop file for the ligand (not necessary for option GC).
    #
    #   GC - Snapshots are generated from trajectories (see below).
    #   AS - Residues are mutated during generation of snapshots from trajectories.
    #   DC - Decompose the free energies into individual contributions
    #        (only works with MM and GB).
    #
    #   MM - Calculation of gas phase energies using sander.
    #   GB - Calculation of desolvation free energies using the GB models in sander
    #        (see below).
    #   PB - Calculation of desolvation free energies using delphi (see below).
    #   MS - Calculation of nonpolar contributions to desolvation using molsurf
    #        (see below).
    #        If MS == 0, nonpolar contributions are calculated with the LCPO method
    #        in sander.
    #   NM - Calculation of entropies with nmode.
    #
    PREFIX                xxx_tr
    PATH                  ./
    #
    COMPLEX               1
    RECEPTOR              1
    LIGAND                1
    #
    COMPT                 compl_tr.prmtop
    RECPT                 prot_tr.prmtop
    LIGPT                 lig.prmtop
    #
    GC                    1
    AS                    0
    DC                    0
    #
    MM                    0
    GB                    0
    PB                    0
    MS                    0
    #
    NM                    1
    #
    ################################################################################
    @MAKECRD
    #
    # The following parameters are passed to make_crd_hg, which extracts snapshots
    #   from trajectory files. (This section is only relevant if GC = 1 OR AS = 1 above.)
    #
    #   BOX - "YES" means that periodic boundary conditions were used during MD
    #         simulation and that box information has been printed in the
    #         trajecotry files; "NO" means opposite.
    #   NTOTAL - Total number of atoms per snapshot printed in the trajectory file
    #            (including water, ions, ...).
    #   NSTART - Start structure extraction from NSTART snapshot.
    #   NSTOP - Stop structure extraction at NSTOP snapshot.
    #   NFREQ - Every NFREQ structure will be extracted from the trajectory.
    #
    #   NUMBER_LIG_GROUPS - Number of subsequent LSTART/LSTOP combinations to
    #                       extract atoms belonging to the ligand.
    #   LSTART - Number of first ligand atom in the trajectory entry.
    #   LSTOP - Number of last ligand atom in the trajectory entry.
    #   NUMBER_REC_GROUPS - Number of subsequent RSTART/RSTOP combinations to
    #                       extract atoms belonging to the receptor.
    #   RSTART - Number of first receptor atom in the trajectory entry.
    #   RSTOP - Number of last receptor atom in the trajectory entry.
    #   Note: If only one molecular species is extracted, use only the receptor
    #         parameters (NUMBER_REC_GROUPS, RSTART, RSTOP).
    #
    BOX                   YES
    NTOTAL                36606
    NSTART                1
    NSTOP                 20
    NFREQ                 1
    #
    NUMBER_LIG_GROUPS     1
    LSTART                2622
    LSTOP                 3862
    NUMBER_REC_GROUPS     1
    RSTART                1
    RSTOP                 2621
    #
    #################################################################################
    @TRAJECTORY
    #
    # Trajectory names
    #
    #   The following trajectories are used to extract snapshots with "make_crd_hg":
    #   Each trajectory name must be preceeded by the TRAJECTORY card.
    #   Subsequent trajectories are considered together; trajectories may be
    #     in ascii as well as in .gz format.
    #   To be able to identify the title line, it must be identical in all files.
    #
    TRAJECTORY            ./mdcrd5
    #
    #################################################################################
    @NM
    #
    # Parameters for sander/nmode calculation (this section is only relevant if NM = 1 above)
    #
    #   The following parameters are passed to sander (for minimization) and nmode
    #     (for entropy calculation using gasphase statistical mechanics).
    #   For further details see documentation.
    #
    #   DIELC - (Distance-dependent) dielectric constant
    #   MAXCYC - Maximum number of cycles of minimization.
    #   DRMS - Convergence criterion for the energy gradient.
    #
    DIELC                 4
    MAXCYC                1000
    DRMS                  0.1
    #
    #################################################################################

    The first mm_pbsa.in1 sample file (pbsa version - does not work with charged ligands - Do not use!)

    It is essentially taken from /home/bio/AMBER/Amber8/src/mm_pbsa/Examples (but files merged and set up to have all files in the same directory).
    You should change issues marked in yellow.

    #
    # Input parameters for mm_pbsa.pl
    # This example just generates snapshots from a trajectory file
    #
    # Holger Gohlke
    # 08.01.2002
    #
    ################################################################################
    @GENERAL
    #
    # General parameters
    #   0: means NO; >0: means YES
    #
    #   mm_pbsa allows to calculate (absolute) free energies for one molecular
    #     species or a free energy difference according to:
    #
    #     Receptor + Ligand = Complex,
    #     DeltaG = G(Complex) - G(Receptor) - G(Ligand).
    #
    #   PREFIX - To the prefix, "{_com, _rec, _lig}.crd.Number" is added during
    #            generation of snapshots as well as during mm_pbsa calculations.
    #   PATH - Specifies the location where to store or get snapshots.
    #
    #   COMPLEX - Set to 1 if free energy difference is calculated.
    #   RECEPTOR - Set to 1 if either (absolute) free energy or free energy
    #              difference are calculated.
    #   LIGAND - Set to 1 if free energy difference is calculated.
    #
    #   COMPT - parmtop file for the complex (not necessary for option GC).
    #   RECPT - parmtop file for the receptor (not necessary for option GC).
    #   LIGPT - parmtop file for the ligand (not necessary for option GC).
    #
    #   GC - Snapshots are generated from trajectories (see below).
    #   AS - Residues are mutated during generation of snapshots from trajectories.
    #   DC - Decompose the free energies into individual contributions
    #        (only works with MM and GB).
    #
    #   MM - Calculation of gas phase energies using sander.
    #   GB - Calculation of desolvation free energies using the GB models in sander
    #        (see below).
    #   PB - Calculation of desolvation free energies using delphi (see below).
    #   MS - Calculation of nonpolar contributions to desolvation using molsurf
    #        (see below).
    #        If MS == 0, nonpolar contributions are calculated with the LCPO method
    #        in sander.
    #   NM - Calculation of entropies with nmode.
    #
    PREFIX                xxx
    PATH                  ./
    #
    COMPLEX               1
    RECEPTOR              1
    LIGAND                1
    #
    COMPT                 compl.prmtop
    RECPT                 prot.prmtop
    LIGPT                 lig.prmtop
    #
    GC                    1
    AS                    0
    DC                    0
    #
    MM                    1
    GB                    1
    PB                    1
    MS                    1
    #
    NM                    0
    #
    ################################################################################
    @MAKECRD
    #
    # The following parameters are passed to make_crd_hg, which extracts snapshots
    #   from trajectory files. (This section is only relevant if GC = 1 OR AS = 1 above.)
    #
    #   BOX - "YES" means that periodic boundary conditions were used during MD
    #         simulation and that box information has been printed in the
    #         trajecotry files; "NO" means opposite.
    #   NTOTAL - Total number of atoms per snapshot printed in the trajectory file
    #            (including water, ions, ...).
    #   NSTART - Start structure extraction from NSTART snapshot.
    #   NSTOP - Stop structure extraction at NSTOP snapshot.
    #   NFREQ - Every NFREQ structure will be extracted from the trajectory.
    #
    #   NUMBER_LIG_GROUPS - Number of subsequent LSTART/LSTOP combinations to
    #                       extract atoms belonging to the ligand.
    #   LSTART - Number of first ligand atom in the trajectory entry.
    #   LSTOP - Number of last ligand atom in the trajectory entry.
    #   NUMBER_REC_GROUPS - Number of subsequent RSTART/RSTOP combinations to
    #                       extract atoms belonging to the receptor.
    #   RSTART - Number of first receptor atom in the trajectory entry.
    #   RSTOP - Number of last receptor atom in the trajectory entry.
    #   Note: If only one molecular species is extracted, use only the receptor
    #         parameters (NUMBER_REC_GROUPS, RSTART, RSTOP).
    #
    BOX                   YES
    NTOTAL                36606
    NSTART                21
    NSTOP                 40
    NFREQ                 1
    #
    NUMBER_LIG_GROUPS     1
    LSTART                2622
    LSTOP                 3862
    NUMBER_REC_GROUPS     1
    RSTART                1
    RSTOP                 2621
    #
    #################################################################################
    @TRAJECTORY
    #
    # Trajectory names
    #
    #   The following trajectories are used to extract snapshots with "make_crd_hg":
    #   Each trajectory name must be preceeded by the TRAJECTORY card.
    #   Subsequent trajectories are considered together; trajectories may be
    #     in ascii as well as in .gz format.
    #   To be able to identify the title line, it must be identical in all files.
    #
    TRAJECTORY            ./mdcrd5
    #
    ################################################################################
    @PB
    #
    # PB parameters (this section is only relevant if PB = 1 above)
    #
    #   The following parameters are passed to the PB solver.
    #   Additional parameters (e.g. SALT) may be added here.
    #   For further details see the delphi and pbsa documentation.
    #
    #   PROC - Determines which method is used for solving the PB equation:
    #          If PROC = 1, the delphi program is applied. If PROC = 2,
    #           the pbsa program of the AMBER suite is used.
    #   REFE - Determines which reference state is taken for PB calc:
    #          If REFE = 0, reaction field energy is calculated with EXDI/INDI.
    #            Here, INDI must agree with DIELC from MM part.
    #          If REFE > 0 && INDI > 1.0, the difference of total energies for
    #            combinations EXDI,INDI and 1.0,INDI is calculated.
    #            The electrostatic contribution is NOT taken from sander here.
    #   INDI - Dielectric constant for the molecule.
    #   EXDI - Dielectric constant for the surrounding solvent.
    #   SCALE - Lattice spacing in no. of grids per Angstrom.
    #   LINIT - No. of iterations with linear PB equation.
    #   PRBRAD - Solvent probe radius in A (e.g. use 1.4 with the PARSE parameter set
    #     and 1.6 with the radii optimized by R. Luo)
    #
    #   Parameters for pbsa only
    #
    #   RADIOPT - Option to set up atomic avity radii for molecular surface calculation
    #     and dielectric assignment. A value of 0 uses the cavity radii from the prmtop file.
    #     A value of 1 sets up optimized cavity radii at the pbsa initialization phase.
    #     The latter radii are optimized for model compounds of proteins only; use cautions
    #     when applying these radii to nucleic acids.
    #
    #   Parameters for delphi only
    #
    #   FOCUS - If FOCUS > 0, subsequent (multiple) PERFIL and SCALE parameters are
    #     used for multiple delphi calculations using the focussing technique.
    #     The # of _focussing_ delphi calculations thus equals the value of FOCUS.
    #   PERFIL - Percentage of the lattice that the largest linear dimension of the
    #            molecule will fill.
    #   CHARGE - Name of the charge file.
    #   SIZE - Name of the size (radii) file.
    #
    #   SURFTEN / SURFOFF - Values used to compute the nonpolar contribution Gnp to
    #                  the desolvation according to Gnp = SURFTEN * SASA + SURFOFF.
    #
    #
    PROC                  2
    REFE                  0
    INDI                  1.0
    EXDI                  80.0
    SCALE                 2.0
    LINIT                 500
    PRBRAD                1.6
    #
    RADIOPT               1
    #
    SURFTEN               0.00542
    SURFOFF               0.92
    #
    ################################################################################
    @MM
    #
    # MM parameters (this section is only relevant if MM = 1 above)
    #
    #   The following parameters are passed to sander.
    #   For further details see the sander documentation.
    #
    #   DIELC - Dielectricity constant for electrostatic interactions.
    #           Note: This is not related to GB calculations.
    #
    DIELC                 1.0
    #
    ################################################################################
    @GB
    #
    # GB parameters (this section is only relevant if GB = 1 above)
    #
    #   The first group of the following parameters are passed to sander.
    #   For further details see the sander documentation.
    #
    #   IGB - Switches between Tsui's GB (1), Onufriev's GB (2, 5).
    #   GBSA - Switches between LCPO (1) and ICOSA (2) method for SASA calc.
    #          Decomposition only works with ICOSA.
    #   SALTCON - Concentration (in M) of 1-1 mobile counterions in solution.
    #   EXTDIEL - Dielectricity constant for the solvent.
    #   INTDIEL - Dielectricity constant for the solute
    #
    #   SURFTEN / SURFOFF - Values used to compute the nonpolar contribution Gnp to
    #                   the desolvation according to Gnp = SURFTEN * SASA + SURFOFF.
    #
    IGB                   2
    GBSA                  1
    SALTCON               0.00
    EXTDIEL               80.0
    INTDIEL               1.0
    #
    SURFTEN               0.00542
    SURFOFF               0.92
    #
    ################################################################################
    @MS
    #
    # Molsurf parameters (this section is only relevant if MS = 1 above)
    #
    #   PROBE - Radius of the probe sphere used to calculate the SAS.
    #           Since Bondi radii are already augmented by 1.4A, PROBE should be 0.0
    #
    PROBE                 0.0
    #
    ################################################################################

Sample mm_pbsa.in3 file for MM/GBSA decomposotion.

#
# Input parameters for mm_pbsa.pl
# This example just generates snapshots from a trajectory file
#
# Holger Gohlke
# 08.01.2002
#
################################################################################
@GENERAL
#
# General parameters
#   0: means NO; >0: means YES
#
#   mm_pbsa allows to calculate (absolute) free energies for one molecular
#     species or a free energy difference according to:
#
#     Receptor + Ligand = Complex,
#     DeltaG = G(Complex) - G(Receptor) - G(Ligand).
#
#   PREFIX - To the prefix, "{_com, _rec, _lig}.crd.Number" is added during
#            generation of snapshots as well as during mm_pbsa calculations.
#   PATH - Specifies the location where to store or get snapshots.
#
#   COMPLEX - Set to 1 if free energy difference is calculated.
#   RECEPTOR - Set to 1 if either (absolute) free energy or free energy
#              difference are calculated.
#   LIGAND - Set to 1 if free energy difference is calculated.
#
#   COMPT - parmtop file for the complex (not necessary for option GC).
#   RECPT - parmtop file for the receptor (not necessary for option GC).
#   LIGPT - parmtop file for the ligand (not necessary for option GC).
#
#   GC - Snapshots are generated from trajectories (see below).
#   AS - Residues are mutated during generation of snapshots from trajectories.
#   DC - Decompose the free energies into individual contributions
#        (only works with MM and GB).
#
#   MM - Calculation of gas phase energies using sander.
#   GB - Calculation of desolvation free energies using the GB models in sander
#        (see below).
#   PB - Calculation of desolvation free energies using delphi (see below).
#   MS - Calculation of nonpolar contributions to desolvation using molsurf
#        (see below).
#        If MS == 0, nonpolar contributions are calculated with the LCPO method
#        in sander.
#   NM - Calculation of entropies with nmode.
#
PREFIX                btn1
PATH                  ./
#
COMPLEX               1
RECEPTOR              1
LIGAND                1
#
COMPT                 compl.prmtop
RECPT                 prot.prmtop
LIGPT                 lig.prmtop
#
GC                    0
AS                    0
DC                    1
#
MM                    1
GB                    1
PB                    0
MS                    0
#
NM                    0
#
#
################################################################################
@MM
#
# MM parameters (this section is only relevant if MM = 1 above)
#
#   The following parameters are passed to sander.
#   For further details see the sander documentation.
#
#   DIELC - Dielectricity constant for electrostatic interactions.
#           Note: This is not related to GB calculations.
#
DIELC                 1.0
#
################################################################################
@GB
#
# GB parameters (this section is only relevant if GB = 1 above)
#
#   The first group of the following parameters are passed to sander.
#   For further details see the sander documentation.
#
#   IGB - Switches between Tsui's GB (1), Onufriev's GB (2, 5).
#   GBSA - Switches between LCPO (1) and ICOSA (2) method for SASA calc.
#          Decomposition only works with ICOSA.
#   SALTCON - Concentration (in M) of 1-1 mobile counterions in solution.
#   EXTDIEL - Dielectricity constant for the solvent.
#   INTDIEL - Dielectricity constant for the solute
#
#   SURFTEN / SURFOFF - Values used to compute the nonpolar contribution Gnp to
#                   the desolvation according to Gnp = SURFTEN * SASA + SURFOFF.
#
IGB                   2
GBSA                  2
SALTCON               0.00
EXTDIEL               80.0
INTDIEL               1.0
#
SURFTEN               0.00542
SURFOFF               0.92
#
################################################################################
@MS
#
# Molsurf parameters (this section is only relevant if MS = 1 above)
#
#   PROBE - Radius of the probe sphere used to calculate the SAS.
#           Since Bondi radii are already augmented by 1.4A, PROBE should be 0.0
#
PROBE                 0.0
#
#
################################################################################
@DECOMP
#
# Energy decomposition parameters (this section is only relevant if DC = 1 above)
#
#   Energy decomposition is performed for gasphase energies, desolvation free
#     energies calculated with GB, and nonpolar contributions to desolvation
#     using the LCPO method.
#   For amino acids, decomposition is also performed with respect to backbone
#     and sidechain atoms.
#
#   DCTYPE - Values of 1 or 2 yield a decomposition on a per-residue basis,
#            values of 3 or 4 yield a decomposition on a pairwise per-residue
#               basis. For the latter, so far the number of pairs must not
#               exceed the number of residues in the molecule considered.
#            Values 1 or 3 add 1-4 interactions to bond contributions.
#            Values 2 or 4 add 1-4 interactions to either electrostatic or vdW
#              contributions.
#
#   COMREC - Residues belonging to the receptor molecule IN THE COMPLEX.
#   COMLIG - Residues belonging to the ligand molecule IN THE COMPLEX.
#   RECRES - Residues in the receptor molecule.
#   LIGRES - Residues in the ligand molecule.
#   {COM,REC,LIG}PRI - Residues considered for output.
#   {REC,LIG}MAP - Residues in the complex which are equivalent to the residues
#                  in the receptor molecule or the ligand molecule.
#
DCTYPE                2
#
COMREC                1-497
COMLIG                498-498
COMPRI                1-498
RECRES                1-497
RECPRI                1-497
RECMAP                1-497
LIGRES                1-1
LIGPRI                1-1
LIGMAP                498-498
#
#
################################################################################


    Errors and solutions

****************************************************************************
 bad atom type: F  (from sander)
From: Holger Gohlke (gohlke_at_scripps.edu)
Date: Wed May 21 2003 - 10:49:20 CDT:
The error occurs because F has not been parameterized for the LCPO model (see Weiser et al., JCC 1999, 20,217), which is used to calculate the solvent-accessible surface area for GBSA calculations within sander. If you "come up" with some reasonable parameters for this, you can add these parameters to the list in mdread.f (in $AMBERHOME/src/sander; see the code after
c --- construct parameters for SA calculation; note that the
c radii stored in L165 are augmented by 1.4 Ang.

I added the following code in Amber9/src/sander/mdread.f (cl actually not needed, but it is treated as carbon otherwise):
! Cl added by UR 5/11-08: radius = 1.8, parameters from S
            else if (atype == 'Cl' .or. atype(1:1) == 'cl') then
               x(l165-1+i) = 1.80d0 + 1.4d0
               x(l170-1+i) = 0.54581d0
               x(l175-1+i) = -0.19477d0
               x(l180-1+i) = -0.0012873d0
               x(l185-1+i) = 0.00029247d0
! F added by UR 5/11-08: Bondi radius and parameters from O nbond=1    
            else if (atype(1:1) == 'F' .or. atype(1:1) == 'f') then
                  x(l165-1+i) = 1.47d0 + 1.4d0
                  x(l170-1+i) = 0.77914d0
                  x(l175-1+i) = -0.25262d0
                  x(l180-1+i) = -0.0016056d0
                  x(l185-1+i) = 0.00035071d0

It turned out that it is also needed for the ICOSA method (gbsa=2; a few lines further down):
! F added by UR 5/11-08 (Bondi); Cl not needed - identical with carbon
            else if (atype(1:1) == 'F' .or. atype(1:1) == 'f') then
               x(L165-1+i) = 1.47d0 + 1.4d0

****************************************************************************
PB Bomb in pb_aaradi(): No radius assigned for atom 4 N1  NT
Added the following code in $AMBERHOME/AmberTools/src/pbsa/sa_driver.f
      if ( radi(iatm) == ZERO ) then
! Changed by UR - if still undefined, use radius from parm file
            radi(iatm) = rin(iatm)
!         write(6, *) 'PB Bomb in pb_aaradi(): No radius assigned for atom', iatm, igraph(iatm), isymbl(iatm)
!         call mexit(6, 1)
      endif

****************************************************************************
PB bomb in pb_reslist(): maxnbr too small
Added the following line (#435) to $AMBERHOME/src/mm_pbsa/mm_pbsa_createinput.pm:

print OUT "  cutres = 12 \n";

****************************************************************************
    PB Bomb in setgrd(): focusing grid too large    2
    reset fillratio to a larger number 2.000

    Added the following line (#435) to $AMBERHOME/src/mm_pbsa/mm_pbsa_createinput.pm:
      print OUT "  fillratio=3 \n";

    ****************************************************************************
    If you have more than 10 NUMBER_REC_GROUPS in the mm_pdbs.in file, you need to change
    C NMO - max number of molecule groups
            parameter (NMO=10)
    in the file $AMBERHOME/src/mm_pbsa/make_crd_hg.f and recompile the program by
    make make_crd_hg

    ****************************************************************************
     PB bomb in pb_setgrd(): Allocation aborted

    You have run out of memory.
    Reduce spacing

    **********************pbsa_com.out_sp******************************************************
    PB bomb in pb_saarc(): Allocation aborted

    You have run out of memory.
    Reduce spacing

    or by using single precision in pbsa and recompile.

    To recompile pbsa in single precision, you can comment out the following
    line in "pb_def.h", i.e. change

    #define PBDPREC

    to

    !#define PBDPREC

    and "make clean" and "make install" again.


    ****************************************************************************
    Missing radii are added in
    mm_pbsa_calceneent.pm



    Molsurf

    Running molsurf as a stand-alone program

    Command:
    molsurf file.pqr radius
    where file.pqr is a pqr file with the coordinates, charges, and radii (pdb file)
    and radius is the solvent probe radius
     According to the mm_pbsa script:
    Bondi radii + 1.4A and probe radius of 0.0A yields SAS
    Bondi radii + 0.0A and probe radius of 1.4A yields molecular surface
    Bondi radii + 0.0A and probe radius of 0.0A yields vdW surface

    Example pqr file:
    ATOM      1 O3   BT1     1       8.139  18.895 -23.268 -0.6  2.8

    The number of decimals is not important, but all fields must be delimited by white-space, even the atom name and residue name.

    pbsa.in file (generated by mm/pbsa)

    File generated by mm_pbsa.pl. Using PB
     &cntrl
      ntf    = 1,       ntb    = 0,
      igb    = 10,      dielc  = 1.0,
      cut    = 999.0,   nsnb   = 99999,
      scnb   = 2.0,     scee   = 1.2,
      imin   = 1,       maxcyc = 0,       ntmin  = 2,
     &end
     &pb
      epsin  = 1.0,     epsout  = 80.0,
      istrng = 0,       radiopt = 1,
      sprob  = 1.6,     space   = 0.5,
      maxitn = 500
      npbverb= 1,
      cutres = 12,fillratio=3
     &end

    sander_com.in file (generated by mm/pbsa)

    The sander_lig.in and sander_rec.in files are identical
    command: sander -i sanmin_com.in -o sanmin_com.out -p compl_tr.prmtop -c btn1_tr_com.crd.1 -r restart

    File generated by mm_pbsa.pl. Using  MM GB
     &cntrl
      ntf    = 1,       ntb    = 0,       dielc  = 1.0,
      idecomp= 0,
      igb    = 2,       saltcon= 0.00,
      offset = 0.09,
      intdiel= 1.0,     extdiel= 80.0,
      gbsa   = 0,       surften= 1.0,
      cut    = 999.0,   nsnb   = 99999,
      scnb   = 2.0,     scee   = 1.2,

      imin   = 1,       maxcyc = 1,       ncyc   = 0,
     &end

    nmode_com.in file (generated by mm/pbsa)

    The nmode_lig.in and nmode_rec.in files are identical
    Command: nmode -i nmode_com.in -o nmode_com.out -p compl_tr.prmtop -c restart

    File generated by mm_pbsa.pl
     &data
      ntx    = 0,
      ntrun  = 1,       nvect  = 0,
      drms   = 1,
      dielc  = 4.0,     idiel  = 0,
      scnb   = 2.0,     scee   = 1.2
     &end
    END

    But if you want to run nmode alone, you normally want to set ntx=1 (input coordinates in a formatted file)