The patch you mentioned was constructed, but its changes were not yet
pushed to the branch as we are awaiting approval. I can easily amend the
patch so that pb_read will no longer force radiopt = 2 when phiout=1 and
will instead display a warning if it is set to some other value.
Be aware that if you use radiopt other than 2, it will be up to you to
find and use the appropriate atomic radii in your visualization software
since the molecular surface will no longer correspond to the one defined
by the radii in your topology file... If you use the levelset description
(current default I believe), you can always print out the level set as a
dx file. If not, you will have to construct the surface yourself.
-Wes
> 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
>
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Dec 03 2013 - 12:30:02 PST
