Re: [AMBER] pbsa: solvation radius from topology file not used when radiopt=0

From: Ray Luo, Ph.D. <ray.luo.uci.edu>
Date: Tue, 3 Dec 2013 10:38:47 -0800

Jan-Philip,

Wes has submitted the patch. I think it's pending approval

Also as for radiopt upon the printout of phi map, we can certainly do it as
you suggested The original argument for hardwiring the radii for
visualization is this. It is not straightforward to set radii in VMD/PyMol
and the phi map is supposed to use values on the outside of the molecular
surface, i.e. solvated phi values since the purpose is to get some
intuition on the effect of a molecular probe approaching the surface from
the outside. Thus we want the best possible match between the radii used in
the visualization software and PBSA. Furthermore. PBSA reduces the radii by
0.1 A to make sure the projected phi values on the molecular surfaces are
solvated values.

However, this is not necessary when the newer DX volumetric data output
from Amber is used for surface rendering. The plan is that we'll still set
radiopt == 2 as the default condition and remind users to reset radii in
the visualization software to those used by Amber if radiopt != 2.

Ray

--
Ray Luo, Ph.D.
Professor,
Biochemistry, Molecular Biophysics, and
Biomedical Engineering
University of California, Irvine, CA 92697-3900
On Tue, Dec 3, 2013 at 4:54 AM, Jan-Philip Gehrcke
<jgehrcke.googlemail.com>wrote:
> Hello Wes and Ray,
>
> thanks for your immediate feedback. I have seen that you have a patch
> since last week -- any news about when this will be made public?
>
> Also, I would like to discuss again the fact that radiopt is internally
> overridden to 2 when phiout == 1. In understand that there were some
> arguments for taking this decision. However, in my opinion, the current
> behavior must be changed for two reasons:
>
> - the current behavior contradicts documentation. The docs say that via
> radiopt = 1 PBSA uses the radii from the prmtop file. When no error
> appears as of a missing atom type, the user has no chance to realize
> that PBSA magically decided to use its hard-coded radii rather than the
> ones from the prmtop.
>
> - one could fix the docs or the code -- I vote for code. There should be
> a way to make PBSA take radii from the prmtop file *without being forced
> to manually change code*. And I don't see why this statement should be
> invalid depending on an output format.
>
>
> Thanks a lot for looking into this,
>
> Jan-Philip
>
>
> On 11/26/2013 07:49 PM, Ray Luo, Ph.D. wrote:
> > Jan-Philip,
> >
> > Okay, it turns out someone overwrote that change in the AmberTools13
> > release. Wes will need to issue a patch to put it back.
> >
> > Thanks a lot for reporting!
> >
> > Ray
> >
> > --
> > Ray Luo, Ph.D.
> > Professor,
> > Biochemistry, Molecular Biophysics, and
> > Biomedical Engineering
> > University of California, Irvine, CA 92697-3900
> >
> >
> > On Mon, Nov 25, 2013 at 10:52 AM, Jan-Philip Gehrcke <
> > jgehrcke.googlemail.com> wrote:
> >
> >> On 11/25/2013 07:13 PM, Ray Luo, Ph.D. wrote:
> >>> Jan-Philip,
> >>>
> >>> Indeed, when you ask for potential printout, pbsa would set radii as
> used
> >>> in VMD for consistency. That's why it overwrites your option in the
> input
> >>> file. If you insist in using radiopt=0, you can certainly comment the
> >> line
> >>> "radiopt=2" output and do it again.
> >>>
> >>> Ray
> >>
> >> Hey Ray,
> >>
> >> thanks for the confirmation. Two things left to be discussed:
> >>
> >> 1) I outcommented radiopt=2. In another attempt I left it as is and
> >> simply added the missing atom types to phi_aaradii(). In both cases,
> >> pbsa runs without error and produces the same result. So far so good.
> >> However, I end up with electrostatic potential values that seem to be
> >> too large by about two orders of magnitude.
> >>
> >> I have used pbsa earlier (1-2 years ago, from AmberTools12), with
> >> basically the same input (radiopt=0 is the only difference here) where I
> >> got values between 1 and 10 kcal/mol for most of the space regions near
> >> the protein surface. Now, I get something between 500 and 2000 (whatever
> >> units). So, the units seem to have changed.
> >>
> >> Back in 2012, we also had an issue with the estatic potential units in
> >> the DX output file: http://archive.ambermd.org/201207/0466.html. At
> >> these times, Wes took care of this and released a bug fix making sure
> >> that the unit would be kcal/(mol*e). Did this accidentally change again?
> >>
> >> I have further proof that something changed. Found the
> >> AmberTools12-based pbsa binary that I used earlier. I used both, the old
> >> binaries and the new binaries with *the same input*. The first few
> >> values of the DX file corresponding to the old binary:
> >>
> >>      1.3067049144E+00    1.3232708113E+00    1.3399177244E+00
> >>      1.3566168815E+00    1.3733524282E+00    1.3900526596E+00
> >>
> >> Using the new AmberTools13-based pbsa binary:
> >>
> >>      5.4732938513E+02    5.5427901231E+02    5.6121241210E+02
> >>      5.6816177213E+02    5.7512228560E+02    5.8210218482E+02
> >>
> >> The factor is about 418.
> >>
> >>
> >> 2) I would suggest to change the current code behavior, because
> >> otherwise what happens is entirely intransparent to the user (we define
> >> radiopt=0 and magically get something undocumented). Would it be
> >> possible to use pymol/vmd-radii only in case radiopt was not explicitly
> >> defined in the input file (as kind of fallback)? Then, one should also
> >> add a corresponding warning/information to the output file, that these
> >> radii will now be applied. A missing-radius-error right after such an
> >> info line would make sense.
> >>
> >>
> >> Thanks for your support,
> >>
> >> Jan-Philip
> >>
> >>
> >>>
> >>> --
> >>> Ray Luo, Ph.D.
> >>> Professor,
> >>> Biochemistry, Molecular Biophysics, and
> >>> Biomedical Engineering
> >>> University of California, Irvine, CA 92697-3900
> >>>
> >>>
> >>> On Mon, Nov 25, 2013 at 10:04 AM, Jan-Philip Gehrcke <
> >>> jgehrcke.googlemail.com> wrote:
> >>>
> >>>> On 11/25/2013 06:38 PM, Jan-Philip Gehrcke wrote:
> >>>>
> >>>>> I used 'set default pbradii bondi' in leap when creating that
> topology
> >>>>> file. I am using AmberTools 13 with the latest patches for pbsa. The
> >>>>> whole pbsa input file is this:
> >>>>>
> >>>>> &cntrl
> >>>>> inp=0,
> >>>>> /
> >>>>> &pb
> >>>>> npbopt=0,
> >>>>> radiopt=0,
> >>>>> npbverb=0,
> >>>>> space=1,
> >>>>> fillratio=1.0,
> >>>>> sasopt=2,
> >>>>> phiout=1,
> >>>>> phiform=2,
> >>>>> npbverb=1,
> >>>>> /
> >>>>>
> >>>>> So, what is the issue here? Why is the given radius not used?
> >>>>>
> >>>>
> >>>> In pb_read.F90 there is this code:
> >>>>
> >>>>       ! set phimap output options when requested
> >>>>
> >>>>       if ( phiout == 1 ) then
> >>>>          outphi = .true.
> >>>>          radiopt = 2
> >>>>          donpsa = .false.
> >>>>          npopt = 0
> >>>>       end if
> >>>>
> >>>>
> >>>> Since I am using phiout=1, radiopt becomes set to the (undocumented)
> >>>> value 2, which in turn makes pb_init() end up in this section:
> >>>>
> >>>>          else if ( radiopt == 2 ) then
> >>>>             call phi_aaradi( natom, isymbl, radi )
> >>>>
> >>>> So, if this quick analysis is right, the question would be why we set
> >>>> radiopt = 2 when phiout == 1.
> >>>>
> >>>>
> >>>> _______________________________________________
> >>>> 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
> >>>
> >>
> >>
> >> _______________________________________________
> >> 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
> >
>
>
> _______________________________________________
> 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 Tue Dec 03 2013 - 11:00:03 PST
Custom Search