Re: [AMBER] Amber 12: CalcError in mmpbsa_py_energy and sander

From: Jan-Philip Gehrcke <jgehrcke.googlemail.com>
Date: Fri, 31 Aug 2012 16:58:31 +0200

On 08/31/2012 04:26 PM, Jason Swails wrote:
> On Fri, Aug 31, 2012 at 9:50 AM, Jan-Philip Gehrcke <jgehrcke.googlemail.com
>> wrote:
>
>> On 08/29/2012 04:17 PM, Jason Swails wrote:
>>> To be honest, I'm not sure where the error is coming from. Error
>> messages
>>> in this case seem to be swallowed up (unless they aren't -- are there any
>>> errors printed out in the _MMPBSA_complex_pb.mdout.0 file?).
>>
>> This is the whole content of this mdout file:
>>
>>
>>> Reading parm file (../../complex_unsolvated.prmtop)
>>> title:
>>>
>>> mm_options: e_debug=3
>>> mm_options: ipb=2
>>> mm_options: inp=2
>>> mm_options: epsin=1.000000
>>> mm_options: epsout=80.000000
>>> mm_options: smoothopt=1
>>> mm_options: istrng=0.000000
>>> mm_options: radiopt=1
>>> mm_options: dprob=1.400000
>>> mm_options: iprob=2.000000
>>> mm_options: npbopt=0
>>> mm_options: solvopt=1
>>> mm_options: accept=0.001000
>>> mm_options: maxitn=1000
>>> mm_options: fillratio=4.000000
>>> mm_options: space=0.500000
>>> mm_options: nfocus=2
>>> mm_options: fscale=8
>>> mm_options: bcopt=5
>>> mm_options: eneopt=2
>>> mm_options: cutnb=0.000000
>>> mm_options: sprob=0.557000
>>> mm_options: cavity_surften=0.037800
>>> mm_options: cavity_offset=-0.569200
>>> Processing ASCII trajectory (_MMPBSA_complex.mdcrd.0)
>>>
>>> Processing frame 1
>>> iter Total bad vdW elect nonpolar EPB
>> frms
>>> eff.c(2660) enb14 --> 481.891
>>> eff.c(2661) eel14 --> 4625.482
>>
>>
>>
>>> There were
>>> significant changes to PBSA, and my guess is you're hitting a bomb in the
>>> code that is forcing a bail-out. Since the PBSA code used by sander and
>>> NAB are the same, I would not expect use_sander to fix this.
>>>
>>> The best bet, IMO, would be to run the energy calc for your complex by
>>> hand, and see if MMPBSA.py is inadvertently squashing the error message
>> (in
>>> which case, this should be fixed). The command for PB should be:
>>>
>>> mmpbsa_py_energy -O -i _MMPBSA_pb.mdin -p <complex_prmtop> -c
>>> _MMPBSA_complex.pdb -y _MMPBSA_complex.mdcrd.0 -o output
>>>
>>> (mdcrd replaced by nc if you use netcdf intermediates).
>>>
>>> I'd be interested to see what you learn.
>>
>>
>> You were right about swallowed error messages. We learned something :-).
>> See:
>>
>>> $ mmpbsa_py_energy -i _MMPBSA_pb.mdin -p ../../complex_unsolvated.prmtop
>> -c _MMPBSA_complex.pdb -y _MMPBSA_complex.mdcrd.0 -o output
>>> PB Bomb in pb_aaradi(): No radius assigned for atom 2021 C1 CG
>>
>
> OK. My intent was to capture printed stderr within Python and redirect it
> back to the system standard error stream.

Good intent :-) Unfortunately, the error was written to stdout.

>
> I've always liked the idea of sending error messages to stderr, although it
> may be helpful to send it to both streams so users know *where* in program
> execution the error occurred (although this should be obvious or
> unnecessary from the error message itself).

You're right. However, I think Amber software is a special case. Due to
the importance of the mdout files, also all error messages should end up
there.

>
> What this error says is that this particular atom, C1 type (CG name, I
> think -- you can use "printDetails .2021" in parmed.py to check the full
> atomic details of this atom), was not assigned a radius.

Output of parmed.py:

> Loaded Amber topology file ../complex_unsolvated.prmtop
>
> Reading input from STDIN...
>> printDetails .2021
>
> The mask .2021 matches 1 atoms:
>
> ATOM RES RESNAME NAME TYPE LJ Radius LJ Depth Mass Charge GB Radius GB Screen
> 2021 126 47Y C1 CG 1.9080 0.1094 12.0100 0.4680 1.7000 0.7200



>
> What force field are you using? Where might this atom type be coming from?
> (e.g., ff99SB, ff99SBildn, ff12SB, ... etc.)
>

This atom type is from Glycam 06 G. From Glycam_06g-1.dat:

> GLYCAM_06_G PARAMETERS (FOR AMBER 8.0, RESP 0.010), COPYRIGHT CCRC 2004
> CG 12.01 ! sp3 C aliphatic


The prmtop file has actually been created using leap from AmberTools 1.5
with FF99SB + custom libs using types from GLYCAM_06_G.

PBSA from Amber 11 / AmberTools 1.5 does not complain about the prmtop
file. How comes that it finds the parameter that PBSA from Amber 12 misses?

> The solution would be to add this atom type explicitly to the PBSA source
> code so it assigns the correct radius. The relevant subroutine is
> pb_aaradi in $AMBERHOME/AmberTools/src/pbsa/sa_driver.F90.

Hm, does this mean that for MMPBSA in Amber 12 all parameters for atoms
from third party force fields need to be hardcoded in that file? I am
sure you don't mean that, right? :)

>
> HTH,
> Jason
>


_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Aug 31 2012 - 08:00:06 PDT
Custom Search