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

From: Ray Luo, Ph.D. <ray.luo.uci.edu>
Date: Mon, 25 Nov 2013 11:16:11 -0800

Jan-Philip,

I remember there was an issue in the unit in the potential printout, it's
supposed to be in kcal/mol-e, but not in the internal unit, and Wes
submitted a bug fix to correct it.

I'm forwarding your email to Wes to take a look.

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
Received on Mon Nov 25 2013 - 11:30:02 PST
Custom Search