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
Received on Tue Dec 03 2013 - 05:00:02 PST
