Thanks so much Bill & Jason for the great help ! The issue has been solved.
Jason, you are correct that MMPBSA.py can stand alone and work well without the help of either receptor_mask or ligand_mask. I greatly appreciate your concise explanations.
Amor
--- On Tue, 24/8/10, Jason Swails <jason.swails.gmail.com> wrote:
From: Jason Swails <jason.swails.gmail.com>
Subject: Re: [AMBER] MMPBSA.py indexError
To: "AMBER Mailing List" <amber.ambermd.org>
Date: Tuesday, 24 August, 2010, 8:30 PM
Hello,
On Tue, Aug 24, 2010 at 7:53 AM, Amor San Juan <amorsanjuan.yahoo.com>wrote:
> I have an input trajectory (PPLW10), such that P=protein,P=peptide,L=ligand
> and W10=10 waters. My task is to do binding energy calculation using
> mmpbsa.py.
>
> Note:
> Protein=1-187
> Peptide=188-199
> Ligand=200
> Water=201-210
>
If you want to see how the Ligand binds to the rest of the complex, then the
receptor_prmtop should have the protein, peptide, and water residues in it.
The complex_prmtop should have all of it, obviously, and the ligand_prmtop
should have just the ligand residue. Otherwise it will not work. Bill is
right -- you need to account for the water in one of your receptor or ligand
prmtop files or set strip_mdcrd=1 and supply a solvated_prmtop and a
complex_prmtop without any waters.
>
> My input script:
> -------
> &general
> verbose=2, interval=1, strip_mdcrd=0,
> receptor_mask=':1-187:188-199', ligand_mask=':200',
>
If your prmtops are designed properly above, you should have no need for
specifying receptor_mask and ligand_mask. MMPBSA.py should be able to
figure it out without any help. You should have MMPBSA.py do this since it
will help catch errors. If it complains about not being able to find the
matching masks, then that tells you something is wrong.
> /
> &gb
> igb=5, saltcon=0.100
> /
>
> EOF
>
> # Execute the program
> $DO_PARALLEL $EXE -O -i mmpbsa.in \
> -cp ../PPLW10.prmtop \
> -rp ../PP.prmtop \
> -lp ../L.prmtop \
> -y ../close.mdcrd \
> -o FIN.out \
>
> -------
>
> Note: L.prmtop contain only the atoms for the ligand, without water. The
> PPLW10.prmtop contains solute (protein+peptide+ligand) and solvent (10
> waters).
>
> After I did MMPBSA.py, below is the error from file.out:
> ----
> Starting gb calculation...
>
> calculating ligand contribution...
> calculating receptor contribution...
> calculating complex contribution...
>
> Calculations complete. Writing output file(s)...
>
> Traceback (most recent call last):
> File "/usr/local/amber10_serial/amber10/bin/MMPBSA.py", line 1397, in
> <module>
> decompout, idecomp, dec_verbose, ligstart, decomprun, surften,
> cavity_surften, temp)
> File "/mnt/home/local/amber10_serial/amber10/src/mmpbsa_py/utils.py",
> line 4595, in PrintFinalResults
> '',finaloutput, debug, numframes, one_trajectory,verbose) # noalascan
> is [avg,stdev] for total DELTA G
> File "/mnt/home/local/amber10_serial/amber10/src/mmpbsa_py/utils.py",
> line 1754, in gboutput
> bonddif[x] = bonddif[x] - bond[x]
> IndexError: list index out of range
>
Differences are tracked by subtracting each term from each frame (complex -
receptor - ligand) for a single trajectory run. If the index is running out
of the range, then this means that there are more complex terms than there
are receptor or ligand terms (so there are fewer receptor or ligand terms,
so eventually "x" will be larger than the total number of terms in the bond
array, hence the error). It may be worth catching this exception and
printing a better error message, but the fact remains that your prmtops are
incompatible.
My best advice would be to take out the receptor_mask and ligand_mask in
your input file and try again.
Hope this helps,
Jason
> ----
>
> And, I checked the _MMPBSA_*mdout files:
>
> MMPBSA_ligand_gb.mdout
> -------
> BOND = 12.8698 ANGLE = 52.2840 DIHED =
> 43.8086
> VDWAALS = -12.0174 EEL = 713.3568 EGB =
> -607.9146
> 1-4 VDW = 8.3908 1-4 EEL = -887.6609 RESTRAINT =
> 0.0000
> ESURF = 4.8120
> minimization completed, ENE= -.67207087E+03 RMS= 0.210204E+02
> -------
> It looks fine to me.
>
> But, the MMPBSA_receptor_gb.mdout
> ---
> BOND = 16559590.1778 ANGLE = 290768.1866 DIHED =
> 10659.1550
> VDWAALS = 11518920.7677 EEL = -4590.3634 EGB =
> -4515.0725
> 1-4 VDW = 3577068.6287 1-4 EEL = 412.9539 RESTRAINT =
> 0.0000
> ESURF = 91.7749
> minimization completed, ENE= 0.31948406E+08 RMS= NaN
> ----
> There is a problem in receptor.
>
> And, finally MMPBSA_complex_gb.mdout
> ----
> BOND = 640.6054 ANGLE = 1594.2719 DIHED =
> 2076.8421
> VDWAALS = -1678.6856 EEL = -14017.8454 EGB =
> -2941.7994
> 1-4 VDW = 767.4366 1-4 EEL = 6808.5420 RESTRAINT =
> 0.0000
> ESURF = 87.9797
> minimization completed, ENE= -.66626527E+04 RMS= 0.181731E+02
> ----
>
>
>
> Is the input file for mmpbsa.py correct ? I know there are ten water
> molecules in the complex (PPLW10.prmtop), and am not sure how to specify it
> so the script can recognize it. Even if there is no problem calculating with
> waters, I can not find the cause why there is index error problem.
>
>
> Amor
>
>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Graduate Student
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Aug 25 2010 - 02:30:03 PDT
