AMBER Archive (2008)

Subject: Re: AMBER: Variable PBCAL values in MM-PBSA calculations

From: Daniel Oehme (oehmes4_at_gmail.com)
Date: Sun Dec 07 2008 - 23:56:58 CST

Ray,

After a weekend of thinking about these issues I think the main
problem is that the pbsa program cannot deal with the extra points of
the TIP5P water model. I believe this because any calculation done
using the EP showed high variance in both the complex and receptor
data, while those without EP, (structural water molecule defined as
TIP3P or deleted) gave satisfactory variances. I am pretty sure that
the only way around this problem is to therefore delete the structural
water molecule and do calculations without it unless anyone else can
suggest otherwise.

In regards to the incorrect definition of radii in the parmtops I dont
think this is a MM-PBSA problem. This is a leap problem as the radii
are defined in leap for when we use radiopt=0 in MM-PBSA calculations.
I originally thought the the radii error only occurred for parmtop
files with only one structural water molecule or for TIP5P waters,
however more thorough analysis of the parmtop which I used to run the
explicit TIP3P simulation shows that the first water molecule has it
hydrogen radii defined incorrectly, with the rest of the radii
defined correctly. Interestingly for the explicit TIP5P water
simulations all hydrogens have their radii defined incorrectly, not
just the first one. Therefore for MM-PBSA calculations using
radiopt=0, and a hybrid approach where at least one explicit water
molecule is defined, there will be a problem.

Thanks,

Daniel Oehme
PhD Student
La Trobe University

On Fri, Dec 5, 2008 at 6:06 AM, Ray Luo <rayhuangluo_at_gmail.com> wrote:
> Daniel,
>
> If you can rule out the issue of radii, the next suspect on the list would
> be that you have very different trajectories for the explicit water between
> your TIP3P and TIP5P runs since you have used different water models in
> PMEMD to begin with. I think it won't be too hard to check this out by
> swapping the prmtop files. Of course the trajectories have to be revised too
> ...
>
>> Why are the radii defined incorrectly in some parmtop's
>> and is there any way of reducing the variability in the PBCAl values
>> in the MM-PBSA calculations?
>
> This is because MMPBSA was originally designed for pure implicit solvent
> which does not require any radii for explicit water molecules. You are
> trying with the so-called hybrid approach, i.e. mix of explicit and implicit
> water, which needs radii for all explicit water molecules. In addition,
> TIP5P is still very new. Most people just started playing with it. So it is
> highly possible the existing code cannot assign radii properly for it.
>
> I remember quite a few Amber users/developers tried the hybrid approach
> before and maybe they can chip in on this.
>
> All the best,
> Ray
>
> ==========================================
> Ray Luo, Ph.D.
> Associate Professor
> Dept Molecular Biology & Biochemistry
> University of California, Irvine, CA 92697
> USPS: PO Box 3900 Email: rluo_at_uci.edu
> Phones: (949) 824-9528, 9562
> Web: http://rayl0.bio.uci.edu/
> ==========================================
>
>
> -----Original Message-----
> From: owner-amber_at_scripps.edu [mailto:owner-amber_at_scripps.edu] On Behalf Of
> Daniel Oehme
> Sent: Wednesday, December 03, 2008 10:35 PM
> To: amber
> Subject: Re: AMBER: Variable PBCAL values in MM-PBSA calculations
>
> Thanks for your feedback Ray.
>
> I was already using radiopt=0 so I tried radiopt=1 but pbsa couldnt
> define radii for the extra points. This could easily be fixed by
> defining radii for EP and recompiling but I thought I should look at
> the radii in the prmtop files that I was using for both the simulation
> and the MM-PBSA calculations first. It turns out that the radii for
> TIP5P in both parmtop files are incorrect with one hydrogen set to
> 0.8, one to 1.2 and the two EP have radii of 1.5 which is the same as
> the oxygen.
>
> Searching through the code I found that radii are defined in unitio.c
> and it is here where a problems lie. The issue with the EP are that
> they are not recognised so they are given the default value of 1.5.
> This therefore effects both the simulation and MM-PBSA prmtop files. I
> dont think this will effect the simulation as the radii are not used
> in sander for PMEMD simulations if I understand things correctly.
>
> However the hydrogen radii issues are different and I cant quite
> figure out why we get different radii. The default hydrogen radii is
> 1.2 but if bonded to an Oxygen it should be 0.8. Therefore it seems
> that only one hydrogen is being recognised as being bonded to oxygen
> with maybe the other one being recognised as being bound to a hydrogen
> as multiply bonded hydrogen used the first bond to define which radii
> to use. The only problem with this theory is that for the TIP3P
> simulation parmtop, the hydrogens are recognized correctly and given
> the correct radii. However for the TIP3P MM-PBSA parmtop, where only
> the structural water is included, the hydrogens are not recognized
> properly and are defined radii as seen for the TIP5P parmtop. It is at
> this point I am confused as to why the TIP3P radii are not recognized
> properly for the structural water parmtop.
>
> Given that I thought I had identified why I was getting highly
> variable PBCAL values, I hand edited the parmtop files to define
> appropriate radii (EP=0.0, hydrogens=0.8) and ran the MM-PBSA
> calculations again. However, I still got the same increased variance
> but this time, the GB energy could not be calculated (value of
> -infinity was given in *all.out files).
> Complex SD Receptor SD
> Ligand SD Delta SD
> PBCAL -2634.95 341.17 -2032.19 66.55
> -45.77 1.23 -556.98 347.61
> GB 0.00 0.00 0.00 0.00
> -39.91 1.26 39.91 1.26
>
> I then tried Ray's suggestion of using the default hydrogen radius of
> 1.3A (even though in unitio.c the default is 1.2).
> Complex SD Receptor SD
> Ligand SD Delta SD
> PBCAL -2206.90 193.60 -2126.48 97.35
> -45.77 1.23 -34.64 216.70
>
> Why are the radii defined incorrectly in some parmtop's
> and is there any way of reducing the variability in the PBCAl values
> in the MM-PBSA calculations?
>
> This is because MMPBSA was originally designed for pure implicit solvent
> which does not require any radii for explicit water molecules. You are
> trying with the so-called hybrid approach, i.e. mix of explicit and implicit
> water, which needs radii for all explicit water molecules. So it is highly
> possible the existing code cannot assign radii properly for water molecules.
>
> Any help with these issues would be greatly appreciated.
>
> Daniel Oehme
> PhD Student
> La Trobe University
>
>
> On Tue, Dec 2, 2008 at 12:55 PM, Ray Luo <rayhuangluo_at_gmail.com> wrote:
>> Daniel,
>>
>> Just my 2c ...
>>
>> If you use radiopt=1, I think the radius optimized for TIP3P may be too
>> small for TIP5P, i.e. the two extra points are too close to the surface of
>> the solvent probe sphere boundary, which may cause highly fluctuating PB
>> energies.
>>
>> Maybe use of radiopt=0 could resolve the problem since the radius for the
>> two hydrogen atoms (1.30A) may be large enough to cover the two EP's.
>>
>> All the best,
>> Ray
>>
>> ==========================================
>> Ray Luo, Ph.D.
>> Associate Professor
>> Dept Molecular Biology & Biochemistry
>> University of California, Irvine, CA 92697
>> USPS: PO Box 3900 Email: rluo_at_uci.edu
>> Phones: (949) 824-9528, 9562
>> Web: http://rayl0.bio.uci.edu/
>> ==========================================
>>
>>
>> -----Original Message-----
>> From: owner-amber_at_scripps.edu [mailto:owner-amber_at_scripps.edu] On Behalf
> Of
>> Daniel Oehme
>> Sent: Monday, December 01, 2008 3:45 PM
>> To: amber
>> Subject: AMBER: Variable PBCAL values in MM-PBSA calculations
>>
>> Hi,
>>
>> I am running some MM_PBSA calculations using AMBER 10 and have been
>> getting some highly variable PBCAL values depending on what water type
>> I assign a structural water molecule in the active site to be. The
>> issue comes when trying to use the TIP5P water type. For TIP3P
>> simulations I get the following results:
>> Complex SD Receptor SD
>> Ligand SD Delta SD
>> PBCAL -2066.21 39.76 -2103.88 39.57
>> -45.30 1.39 82.96 7.24
>>
>> while for the TIP5P I get:
>> Complex SD Receptor SD
>> Ligand SD Delta SD
>> PBCAL -2549.08 369.70 -2188.62 295.53
>> -45.77 1.23 314.69 243.66
>>
>> Please note that two different simulations were run to get these
>> results. As can be seen, the Complex and Receptor have extremely high
>> SD's but the Ligands values are fairly similar. If however I run the
>> MM-PBSA calculations without the structural water molecules in the
>> TIP5P simulation I get the following results:
>> Complex SD Receptor SD
>> Ligand SD Delta SD
>> PBCAL -1982.99 66.10 -2024.60 66.13
>> -45.77 1.23 87.38 4.61
>>
>> I have also run a calculation on the TIP5P simulation by making the
>> structural water molecule TIP3P instead of TIP5P and got these
>> results:
>> Complex SD Receptor SD
>> Ligand SD Delta SD
>> PBCAL -1982.75 66.07 -2027.20 65.99
>> -45.77 1.23 90.23 4.65
>>
>> As can be seen, the last two sets of results give SD's are a lot lower
>> in this case yet still a bit higher than in the inital TIP3P
>> simulations. My theory is that the PBCAL calculations are getting
>> confused by the extra points (EP) that are included in the TIP5P water
>> model compared to TIP3P. My question is can anybody think of a reason
>> as to why the PBCAL values would get confused with the EP's and how I
>> can change the setup of the calculation to be able to run the MM-PBSA
>> calculations with a TIP5P structural water molecule?
>>
>> MM-PBSA calculations were run without a solvent sphere and with
>> standard input except for fillratio which was increased to 4.
>>
>> Daniel Oehme
>> PhD Student
>> La Trobe University
>> -----------------------------------------------------------------------
>> The AMBER Mail Reflector
>> To post, send mail to amber_at_scripps.edu
>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>> to majordomo_at_scripps.edu
>>
>> -----------------------------------------------------------------------
>> The AMBER Mail Reflector
>> To post, send mail to amber_at_scripps.edu
>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>> to majordomo_at_scripps.edu
>>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
> to majordomo_at_scripps.edu
>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
> to majordomo_at_scripps.edu
>
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
      to majordomo_at_scripps.edu