Custom Search
-- 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/amberReceived on Tue Nov 26 2013 - 11:00:02 PST