Re: Free energy calculation yields strange results

From: Jun Zhang (coolrainbow_at_yahoo.cn)
Date: Tue May 11 2010 - 21:25:57 CDT

Hi Jérôme Hénin:

Thank you very much for your reply! I have tried the FEP in vacuo and I do get a ZERO result by the following arguments. Does this mean there is not any problems in my force field parameters?

alchVdwLambdaEnd 1.0
alchElecLambdaStart 0.5
alchVdWShiftCoeff 6.0
alchDecouple yes
alchEquilSteps 5000
run 10000

My target molecule in R(aq)--S(aq) FEP calculation involves 16 atoms disappearing and 16 atoms appearing. To test my configuration, I had simplified the molecule(replace a long chain by a hydrogen) and had calculated the FEP of the process

miniR(aq) -- miniS(aq)

, which involves only 3 atoms disappearing and 3 atoms appearing, in aqueous. What surprised me is that without a soft-core potential I can get a ZEOR result:

WITH soft-core potential:

lchVdwLambdaEnd 1.0
alchElecLambdaStart 0.5
alchVdWShiftCoeff 6.0
alchDecouple yes
alchEquilSteps 1000
set numSteps 1000

window net change
0.0 0.0
0.05 -0.4193
0.1 -0.2837
0.15 -0.0038
0.2 0.3774
0.25 0.6252
0.3 1.4203
0.35 1.2897
0.4 1.3728
0.45 1.1525
0.5 1.1773
0.55 1.0678
0.6 1.0874
0.65 1.4239
0.7 1.6932
0.75 1.9875
0.8 2.1508
0.85 2.5398
0.9 3.2771
0.95 3.538
1 3.8044

WITHOUT soft-core potential:

alchVdwLambdaEnd 1.0
alchElecLambdaStart 1.0
alchVdWShiftCoeff 0.0
alchDecouple no
alchEquilSteps 1000
set numSteps 1000

window net change
0.0 0.0
0.05 0.4786
0.1 0.5852
0.15 0.7195
0.2 0.8346
0.25 0.8222
0.3 0.7789
0.35 0.792
0.4 0.9734
0.45 1.0781
0.5 1.0515
0.55 1.1112
0.6 0.9862
0.65 1.0949
0.7 1.0057
0.75 0.8328
0.8 0.6884
0.85 0.5455
0.9 0.5584
0.95 0.0632
1 -0.1426

This is strange because in my opinion soft-core potential usually yields better results.

Then I tried my target molecule (involving 16 appearing and 16 disappearing atoms), however, neither with nor without soft-core potential works:

window with soft without soft
0.0 0.0 0.0
0.05 68.0538 66.6033
0.1 138.18 125.502
0.15 209.085 183.996
0.2 279.447 242.828
0.25 349.937 301.042
0.3 418.634 359.287
0.35 485.52 417.705
0.4 550.931 476.015
0.45 613.22 534.227
0.5 672.778 592.173
0.55 730.448 649.799
0.6 786.848 707.364
0.65 841.011 765.282
0.7 893.274 823.226
0.75 943.493 880.766
0.8 991.912 938.088
0.85 1040.59 995.253
0.9 1090.46 1052.67
0.95 1139.29 1109.44
1 1187.51 1165.92

Increasing alchVdWShiftCoeff does not work. Is the 16 atoms change too steep for the alchemical transformation so causes a terrible sampling problem? I am performing FEP simulation with a longer timescale and waiting to see the results.

Thank you for your suggestions!

Jun Zhang
Nankai University
coolrainbow_at_yahoo.cn

--- 10å¹´5月10日，周一, Jérôme Hénin <jhenin_at_ifr88.cnrs-mrs.fr> 写é“：

> å‘件人: Jérôme Hénin <jhenin_at_ifr88.cnrs-mrs.fr>
> 主题: Re: namd-l: Free energy calculation yields strange results
> 收件人: "Jun Zhang" <coolrainbow_at_yahoo.cn>
> 抄é€: "NAMD MAILLIST" <namd-l_at_ks.uiuc.edu>
> 日期: 2010å¹´5月10æ—¥,周一,下åˆ9:04
> Hi,
>
> The most obvious problem here is the energy minimization
> before each
> FEP stage. It is not recommended, as 500 MD steps will not
> be
> sufficient to thermalize the system back to room
> temperature. There
> should be no need for minimization there.
>
> In addition, there could possibly be a mistake in your
> setup, or a
> (large) convergence/sampling problem. Interactions between
> ions are
> notoriously difficult to sample:
> http://www.ncbi.nlm.nih.gov/pubmed/15584080
>
> That is all I can say based on the information you
> provided. One test
> you can run is the same transformation in vacuum (remove
> the PBCs and
> forget about neutralizing the charge).
>
> Cheers,
> Jerome
>
> On 9 May 2010 14:53, Jun Zhang <coolrainbow_at_yahoo.cn>
> wrote:
> > Hi Everyone!
> >
> > I am perfoming alchemical free energy calculation
> using NAMD 2.7b2. The transformation of the species is from
> chirality R to S
> >
> > A(R)  ---   A(S)
> >
> > in TIP3P water molecules. The box is about 50A. The
> species has a COO and I have use some ions (Na+ and Cl-) to
> neutralize it. I think the free energy should be ZERO since
> water is achiral. However, my result is 2328.95
> > kcal/mol, which made me surprised. My alchemical
> parameters set is:
> >
> > alch             
>       on
> > alchType           
>     FEP
> > alchFile           
>     R2S_wb.fep
> > alchCol           
>      B
> > alchOutFile         
>    R2S_forward.fepout
> > alchOutFreq         
>    10
> >
> > alchVdwLambdaEnd        1.0
> > alchElecLambdaStart     0.5
> > alchVdWShiftCoeff       6.0
> > alchDecouple           
> yes
> >
> > alchEquilSteps          500
> > set numSteps           
> 5000
> >
> > set Lambda0 0.0
> > set dLambda 0.05
> >
> > while {$Lambda0 <= 1.0} {
> >
> > alchLambda $Lambda0
> > set Lambda0 [expr \$Lambda0 + \$dLambda]
> > alchLambda2 $Lambda0
> > minimize 500
> > run $numSteps
> >
> > }
> >
> > Can anyone give some suggestions? Thank you in
> advance.
> >
> >
> > Jun Zhang
> > Nankai University
> > coolrainbow_at_yahoo.cn
> >
> >
> >
> >
> >
>

      

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:54:07 CST