注册 登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

计算生物学实验室

http://www.bioms.net

 
 
 

日志

 
 

Gibbs (free energy) questions  

2010-12-08 19:10:56|  分类: AMBER |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

Gibbs setup issues


Leap-generated input


Defining states

My biggest concern is to the input of the two states that will be compared by the program. The literature states that the 0 state is the coordinates given by PARM, and the 1 state is contained within the PREP file. Yet the description of the command line for GIBBS has no input PREP file.

The word 'state' refers to the _interpretation_ given to the coordinates, i.e. the nature of the atoms, bonds, etc. There is only one set of coordinates, which is normally a restrt file from sander or gibbs. The PREP 'state' is carried via LINK and EDIT to the PARM topology file (in leap, it goes directly into the topology, and the perturbed state is defined via the xleap 'unit Editor'). From the gibbs manual:

       The  state  defined  in PREP is the =1 state and the state    given in the PARM is the =0  state.   The  default  state  from    which  to  start  the  perturbation  is usually =1, because the    coordinates which are carried from EDIT to PARM to MIN to GIBBS    will correspond to the PREP state.    

Can I perturb by having two copies of a residue that don't "see" each other, disappearing one while appearing the other?

Yes; it's complicated though..


Shrinking distances w/ dummy atoms

I am perturbing an ethyl group to nothing. The dummy atoms are separated by a min distance of 0.001. The hydrogens are too close to the other hydrogens and there are bad contacts. However I have at frcmod.dat the nonbonded parameters very low for the atoms ...
              R*    epsilon           DC  0.001   0.0           DH  0.001   0.0  
so even though they are very close are much smaller in size but they still create bad contacts.

I set the dummy atoms closer than the corresponding real atoms so when they disappear there will not be a big "cap" in the system.

The problem here is that MD trajectories almost always become very unstable if you try to shrink the distances to dummy atoms to extremely small distances. This has to do with increasing vibrational frequencies of short length bonds. In my experience, 0.2 A is typically around the shortest distance you can shrink a bond to without encountering problems.

If you simply insist on shrinking to a shorter distance, you will need to reduce the timestep--possibly to something considerably smaller than 0.5 fsec...Given the tremendously decreased trajectory efficiency that results (i.e. you have to run a lot long to perform the same number of psec of dynamics), you are usually better off with the (possibly) somewhat diminished sampling efficiency associated with shrinking to a larger bond length like 0.2 A.


Interpreting results


Sign of energy; forward/back

> Given the following:  >   >  (1) I have defined molecule X in the PREP input.  >  (2) I have defined molecule Y in the PARM input.  >  (3) I have equilibrated two separate systems based on this same  >      set of input files.  One was equilibrated at Lambda=1 (X) and the   >      other was equilibrated at Lamdba=0 (Y).  >   > If I then perform a Fixed-width window Free Energy Perturbation on each  > equilibrated system separately (one X[lam=1] -> Y[lam=0] perturbation,   > one Y[lam=0] -> X[lam=1] ), I get "forward" free energy differences which  > are *both* negative.  >   > Is it true (as it seems to be based on the info on pg 107 of the manual)   > that the *sign* of both energies will correspond to that of the Y -> X   > (lam=0 -> lam=1) perturbation?  In other words, is it true that *real*   > value of delta G for the X -> Y (1->0) is actually the *negative* of the   > reported "forward" value found in the output?  
Yes. You are correct. The "FORWARD" free energy reported is ALWAYS the value corresponding to the change (0->1), regardless of the actual direction in which the simulation was performed. So if you want the free energy for the (1->0) change, take the negative of the energy reported. (NOTE that this applies to Amber 4.0 and later versions. In earlier (3A, 3.0) versions of Gibbs, the meaning of "Forward" for FEP depended on the direction in which the simulation was run.)

Single vs. double-wide trajectories

I decided to try to save computer time by setting the IDWIDE parameter to 1 in order to turn off the double wide sampling ... I ran a test gas phase perturbation. ... the reported final energy differences were not the same for the two calculations. ... Why should these values differ?

Setting IDWIDE = 1 changes the trajectory, and hence the computed free energy (although values calculated from any trajectory will be the same in the limit of infinite sampling, in the finite sampling domain, they will differ). To understand, when IDWIDE=0 (double wide sampling), the trajectory is calculated at the following points:

        0  dl   2dl   3dl   4dl   ...   1    and       dG (forward) = dG (0->dl) + dG (dl->2dl) + dG (2dl->3dl) + ...      dG (reverse) = dG (dl->0) + dG (2dl->dl) + dG (3dl->2dl) + ...    
whereas if IDWIDE = 1, then the trajectory is calculated at the following points:
           dl        3dl        5dl      ...    and       dG (forward) =               -dG (dl->0) + dG (dl->2dl) - dG (3dl->2dl) + dG (3dl->4dl) + ...      (You don't get a "reverse" free energy).    
Note that since are running MD at a different set of "lambda" points, the trajectory will be different, and the net free energy will be different. The difference will be small, of course, if the free energy is converged.

Dihedral angle warnings

 What might be the reason for the warning that the dihedral    angles are not between 0 or 180?               nb-update: total= 186068; regular= 96590; h-b= 77272; pert= 12206            Warning (PEPHI): phase angle not 0 or 180 for torsion:               0- 3- 12- 15 N= 2.0000 Phi= -1.9722            Warning (PEPHI): phase angle not 0 or 180 for torsion:               0- 3- 12- 42 N= 2.0000 Phi= 1.1519    This message is telling you that some of your torsional parameters have  been defined with a phase angle other than the standard values of 0 or 180.  That's phi in the equation:        Ephi = K (1 + cos (N tau - phi)).     The only reason you'd typically use a phi other than 0 or 180 is to  impose a particular torsional restraint. Maybe you're doing this.   (However, I'd suggest that if you are, that you shouldn't really be  doing it in the parm file, but rather with added restraints in the  Sander input). Anyway...    The integers refer to the atom numbers of the affected torsion. However,  you need to convert them to the actual atom numbers using the formula        atom-number = (I/3) + 1    By figuring out what the atom type parameters are for this torsion, you  should be able to quickly locate the appropriate torisonal parameter in the  parameters file. Be sure to observe that the dihedral parameter input  is formatted, with a format:        IPT , JPT , KPT , LPT , IDIVF , PK , PHASE , PN        FORMAT(A2,1X,A2,1X,A2,1X,A2,I4,3F15.2)    (as described in the PARM manual).    If you didn't intend for the phase angles given above, check to see that  your input falls in the correct columns...    

Unexpected delta G = 0

I am trying to calculate the free energy difference of CHCl3 mutated to CCl4 inside some cryptates. I just tried the mutation of CHCl3 ALONE to CCl4 in vacuo. Why is the free enrgy I get identically zero?

In CHCl3 (or ChCl4) you have no non-bonded interactions (remember no non-bonds for 1-3 pairs), so there will be no non-bonded contribution. There are no dihedrals. So the only contributions will be from bond and angle changes. To get these calculated correctly and reported, you need to set INTPRT=5 (report intra-pert-group contributions), NTC=NTF=3 (turn SHAKE on for the solute) and ICORC=1 (calculate the "PMF Bond contribution").


Hysteresis & center of mass motion

In perturbing from ACE-GLY-NME dipeptide to ACE-ALA-NME in gas phase, the free energy values show a large hysteresis between the forward and backward runs (double-wide sampling). On the other hand, if I start with the alanine dipeptide as the initial residue and perturb to the glycine dipeptide, there is little hysteresis.

This large hysteresis for the former run disappears if I simply do *not* remove the translational and rotational motion from the system. This problem does not seem to arise if the perturbations are done in a solvated environment, or in a full protein environment (i.e., GLY and ALA are part of a folded protein).

As I have noted in a few papers, one should never remove COM motion when performing free energy calculations. This will definitely have adverse affects on the calculated free energies. I belive this arises because the moment of inertia of the molecule is changing, which will affect the partitioning among particular translational/rotational/vibrational modes. Removing COM motion messes this up.


PMF: the CORC term

> I'm doing a PMF calculation moving an atom in between two dummy  > atoms in the way that the distance atom-dummy1 increases while  > the distance dummy2 decreases. In order to do so I have to use  > the option NCORC=1. I now would like to understand the meaning of  > the CORC value printed in the listing file and its relation to  > the free energy value.  
The "CORC" value is, in principal, the contribution from the constrained bonds to the total free energy. With FEP, the components are only approximate in the best of cases, and for various reasons, a lot of non-bonded interactions get lumped into "CORC" in the output. If you're running FEP, you should definitely ignore the CORC value and focus on the total free energy.

With TI, the CORC contribution is better deconvoluted from the rest. But once again, it is not really a good idea to focus on this component alone for most purposes. I'd suggest you still focus primarily on the total free energy.


"Coordinate coupling" = PMF

How can we carry out coordinate coupling with the program GIBBS on AMBER 41 ? This coupling was mentioned in a publication : Rao, B. and al. JACS, vol 111, 9, 1989, 2135.

The procedure that is termed "Coordinate Coupling" in that paper corresponds to the "PMF Bond Contribution," which was detailed in the paper: Pearlman & Kollman (1991) JCP 94, 4532-4545.

The analogue of this contribution, but for TI calculations, is described in Pearlman (1994) JPC 98, 1487-1493.

To include this contribution in your calculations, you set NCORC=1 and NTC=NTF=3 in the Gibbs input. The program takes care of the rest for you automatically.

  评论这张
 
阅读(1370)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018