Re: [AMBER] MMPBSA.py indexError

From: Amor San Juan <amorsanjuan.yahoo.com>
Date: Wed, 25 Aug 2010 02:09:53 -0700 (PDT)

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
Custom Search