Custom Search
-- 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/amberReceived on Tue Dec 03 2013 - 11:00:03 PST