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