Hello,
Can you send the files you are running pbsa on (inpcrd prmtop and pbsa.in
files) so I can take a look and replicate the results. This should help
me figure out what is going wrong.
-Wes
> 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
>
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Nov 25 2013 - 12:30:02 PST
