AMBER Archive (2005)

Subject: Re: AMBER: TI-FEP for VAL --> ALA

From: Ilyas Yildirim (yildirim_at_pas.rochester.edu)
Date: Sat Jun 25 2005 - 18:39:11 CDT

Hey Jiten,

On Sun, 26 Jun 2005, Jiten wrote:

> Hi Ilyas Yildirim,
>
> Thanks for the reply. First of all I apologize, what I have got zero
> difference is of V--> A TI-FEP. The frcmod.DH was constructed as give before
> for ALA-->GLY.

I assume by TI-FEP u mean free energy calculation using Thermodynamic
Integration approach.

>
> First - I do not understand why I need use hydrogen mass for dummy atom.
> Second - I think that if I put the force constant and eq distance/angle for
> the dummy atoms, with that just replaced by that off
> !CT-HC, !HC-CT-HC,!HC-CT-HC etc -- I am confused if it still be alright?
>

In the TI Approach tutorial, the dummy atoms had some mass. I think the
reason is because of the mixed (hybrid) Hamiltonian. The forces are
calculated according to this Hamiltonian. A zero mass might give an error,
I think. The best choice is to use the mass of DH as

MASS
DH 1.008

If it was a transformation from C->DC (dummy Carbon), it is better to use
the mass of carbon for the mass of DC.

For the second part of the question: I would do the same; defining the
missing parameters by making some analogies with the old structure. In
your case you have to define one more parameter, which is the dihedral
angle of

X -CT-HO-X

Now, the reason why you have to have these force constants is to keep the
structure in shape. The question is: Are these extra dummy atom parameters
going to change the energy of the system? Yes they will but I do not think
this is that crucial. These parameters are going to give some extra
internal energy to the system, namely energy coming
because of bonded parameters. In order to compare this free energy change
with an experimental data, you have to do another free energy change
simulation in gas phase (assuming you have done the first simulation in
water). And then, u will subtract these 2 results. These internal/bonded
energies will cancel at the end. So, that result will give you some
reasonable prediction which can be compared with an experimental data.

Here is what I mean:

 C-(H)3 ------------------> H-(DH)3
(in water) G3 (in water)
   /\ /\
   || ||
   || ||
   || ||
   || G1 G2 ||
   || ||
   || ||
   || G4 ||
 C-(H)3 ------------------> H-(DH)3
(in gas) (in gas)

G1 and G2 are the experimental results, which are the solvation free
energies. G3 and G4 are supposed to be the free energy simulation results.
So,

G1 - G2 = G3 - G4

What I mean is, you have to do two free energy calculation, one for the
solvated case, one in gas phase in order to compare it with an
experimental result.

> Also, I expect that the total energy after 12 windows be lower than the
> first window (similar to that of MM-PBSA), but they are not as given below.
> The structure after 12th window is still alright - nothing unusual.
> Therefore, I feel that there something very wrong with it.
>
> Amber TI-FEP expert would help me ???
>
> MM-PBSA A V
> ELE -8.44 -8.02
> VDW -36.38 -36.02
> GAS -44.82 -44.04
> PBSUR -3.14 -3.01
> PBCAL 23.04 24.59
> PBSOL 19.91 21.58
> PBELE 14.61 16.58
> PBTOT -24.91 -22.46
> TSTOT 18.48 17.58
> PBTOT-S -6.43 -4.88
>
>
> This is in good agreement with experiment.
>
> From first window,
>
> A V E R A G E S O V E R 50000 S T E P S
>
>
> NSTEP = 100000 TIME(PS) = 120.000 TEMP(K) = 299.17 PRESS =
> 0.0
> Etot = -77963.7518 EKtot = 19966.8759 EPtot
> = -97930.6277
> BOND = 722.3367 ANGLE = 1873.6537 DIHED =
> 2212.5509
> 1-4 NB = 935.9809 1-4 EEL = 11755.7341 VDWAALS =
> 11246.8856
> EELEC = -126677.7697 EHBOND = 0.0000 RESTRAINT =
> 0.0000
> DV/DL = 768.0586
> Ewald error estimate: 0.6501E-04
> ------------------------------------------------------------------------------From 12th window, A V E R A G E S O V E R 50000 S T E P S NSTEP = 100000 TIME(PS) = 1440.000 TEMP(K) = 299.59 PRESS =0.0 Etot = -77806.2969 EKtot = 19994.8561 EPtot = -97801.1530 BOND = 712.4664 ANGLE = 1866.1103 DIHED =2207.5175 1-4 NB = 933.9198 1-4 EEL = 11775.2910 VDWAALS =11186.6417 EELEC = -126483.0998 EHBOND = 0.0000 RESTRAINT =0.0000 DV/DL = 0.0000 Ewald error estimate: 0.6711E-04 ------------------------------------------------------------------------------The dV/dl are as follows, 1 2 3 4 5 6 7 8 9 10 11 12 lam(i) 0.00922 0.04794 0.11505 0.20634 0.31608 0.43738 0.562620.68392 0.79366 0.88495 0.95206 0.99078 w(i) 0.02359 0.05347 0.08004 0.10158 0.11675 0.12457 0.12457 0.116750.10158 0.08004 0.05347 0.02359 768.05860 361.0522012.45120 -118.71170 -113.01570 -65.58740 -2
 8!
> .53030 -10.29490 -2.07680 -0.17480 -0.00800 0.00000 Total 18.11850 19.305460.99659 -12.05873 -13.19458 -8.17022 -3.55402 -1.20193 -0.21096 -0.01399 -0.00043 0.00000 0.01569Sincerely,Jiten----- Original Message -----From: "Ilyas Yildirim" <yildirim_at_pas.rochester.edu>To: <amber_at_scripps.edu>Sent: Saturday, June 25, 2005 3:48 PMSubject: Re: AMBER: TI-FEP for ALA --> GLY> Hi Jiten,>> The frcmod file u created looks ok, but why do u have 0 mass on the dummy> atoms? I would have chosen the hydrogen mass rather than a zero mass for> the DH atom.>> Second thing, did u check out one of the trajectory file to see if> everything is ok? And also, what are your dv/dl results for these 12> lambda values. I have done a simple free energy calculation for> ethane->methane which looks similar to what you have done (-CH3 changes to> H-DH3). I also got a result which was close to zero. I did the simulation> just for a test purpose. I have used explicit solvent. In order to compare> with s
 o!
> mething real (with an experimental data), u have to do 2 different> TI
> simulation: One in a solvent, the other in gas phase, and take the> difference of these results. That will give u a prediction for an> observable. (Thermodynamic Cycle)>> For the ethane->methane case, the solvation of ethane and methane are> almost the same, so it is reasonable to expect something close to zero.> What is the solvation free energies of ALA and GLY (or the solvation free> energies of the initial and final states of your system)? The difference> of these energies should be equal to the free energy change of the> transformation in a solvent subtracted from the free energy change of the> transformation in gas phase.>> Good luck,>> On Sat, 25 Jun 2005, Jiten wrote:>>> Dear Amber users>>>> For my protein-drug system : while doing the free energy peturbationusing TI method using 12 windows, and using the perl script attached here, Igot almost zero difference in free energy. I used the following atom types,where DH is dummy atom with zero charge and the charges of
 t!
> he alaline isreplaced by glycine charges in the pert.charge options>>>>>> H1 HC H0 DH>> | | | |>> --CT---CT--HC ---> --CT --H0--DH>> | |>> HC DH>>>> My frcmod.DH is as follows,>>>> ALA to GLY>> MASS>> DH 0.000>>>> BOND>> H0-DH 340.0 1.090 !CT-HC>>>> ANGLE>> DH - H0 - DH 35.0 109.50 !HC-CT-HC>> CT - H0 - DH 35.0 109.50 !HC-CT-HC>>>> DIHE>> - The default values as given by leap>> -->> -->>>> NONB>> DH 0.000 0.00>>>> Am I doing something stupid here ? Do I need to set the Bond and angles :K and R values as zero ?>>>> However, I calculate the absolute free energy of the two syatems usingMM-PBSA from the last 1ns trajectory of the 2ns MD simulations (wellconverged - rmsd, temp and energy) and I got that there is clear differencesbetween the two protein systems.
>!
> >>> ELE -8.99 -8.96>> VDW -41.125 -37.68>> GAS
> -50.115 -46.64>> PBSUR -3.025 -3.07>> PBCAL 27.105 24.16>> PBSOL 24.09 21.09>> PBELE 18.12 15.2>> PBTOT -26.025 -25.55>> TSTOT 12.675 15.79>> PBTOT-S -13.35 -9.76>>>>>> Thanks for any suggestions,>>>> Sincerely,>>>> N. Jiten Singh>> C/O Prof. Kwang S. Kim>> Department of Chemistry>> Pohang University of Science and Technology>> San 31, Hyojadong, Namgu>> Pohang 790-784, Korea>> Phone : 82-54-279-5853 ( Lab ) / 279-4138 ( Appt )>> Fax : 82-54-279-8137 (or +82-54-279-3399)>> Web : http://csm.postech.ac.kr/ and http://www.postech.ac.kr/e>> Home Page : http://www.geocities.com/njs_19>> --> Ilyas Yildirim> ---------------------------------------------------------------> - Department of Chemisty - -> - University of Rochester - -> - Hutchison Hall, # B10 - -> - Rochester, NY 14627-0216 - Ph.:(585) 275 67 66 (Office) -> - http://www.pas.rochester.edu/~yildirim/ -> ------------------------------------
 -!
> -------------------------->> -----------------------------------------------------------------------> The AMBER Mail Reflector> To post, send mail to amber_at_scripps.edu> To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu>>>
>

I wrote a small script to calculate the free energy change of this
simulation, and yes, it gave a deltaG = 0.15 (approx). But as I said, this
result mean nothing. U have to do another simulation and subtract it with
this one.

If the only change is -CH3 --> -H, maybe it is reasonable. The
experimental free energy of hydration of methane is 2.0 kcal/mol while it
is 1.8 kcal/mol for ethane. The difference is 0.2 kcal/mol. (JACKS
1996, 118, 6285-6294; original data from Cabani et al. 1981) So, maybe the
same thing is true for your system. Just a thought.

Good luck,

>
>
>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu
>
> 

-- 
  Ilyas Yildirim
  ---------------------------------------------------------------
  - Department of Chemisty       -				-
  - University of Rochester      -				-
  - Hutchison Hall, # B10        -				-
  - Rochester, NY 14627-0216     - Ph.:(585) 275 67 66 (Office)	-
  - http://www.pas.rochester.edu/~yildirim/			-
  ---------------------------------------------------------------

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu