Re: FEP yields the same results for forward and backward computation!

From: Jun Zhang (coolrainbow_at_yahoo.cn)
Date: Thu May 13 2010 - 08:48:42 CDT

Hello:

Thank you for your reply. I have just studied FEP for several days but I have encountered many problems. The "minimization", although I know it should appear, appears because without it I will get an infinite energy:

FEP: RESETTING FOR NEW FEP WINDOW LAMBDA SET TO 0 LAMBDA2 0.1
188 FEP: WINDOW TO HAVE 100 STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.
FEP: USING CONSTANT TEMPERATURE OF 298 K FOR FEP CALCULATION
PRESSURE: 0 nan nan nan nan nan nan nan nan nan
GPRESSURE: 0 nan nan nan nan nan nan nan nan nan

however, with one minimization before each window will eliminate this problem. I don't understand why. And in fact, I have tried very long run but I still can't get a reasonable result.

In my opinion, even if I use a short run(my system is very small), the results of the forward and backward should at least have an opposite sign, but my calculation yielded the same sign. Maybe there exists some stupid problems in my configuration file or system setup. I am still looking for it. Thank you for your ideas.

Yours sincerely

Jun Zhang
Nankai University
coolrainbow_at_yahoo.cn

--- 10年5月13日,周四, Felipe Merino <felmerino_at_uchile.cl> 写道:

> 发件人: Felipe Merino <felmerino_at_uchile.cl>
> 主题: Re: namd-l: FEP yields the same results for forward and backward computation!
> 收件人: "Jun Zhang" <coolrainbow_at_yahoo.cn>
> 抄送: "NAMD MAILLIST" <namd-l_at_ks.uiuc.edu>
> 日期: 2010年5月13日,周四,下午2:24
> Hey,
>
> What you said in the email does not make much sense, but i
> suppose that
> your concern is about the difference between the forward
> and backward
> directions.
>
> I don't know why you make a minimization before every
> window.
>
> Besides that you are using an extremely short calculation
> time per
> window which of course could explain the huge hysteresis
>
> best
>
> felipe
>
> Jun Zhang escribi贸:
> > Hi Everyone!
> >
> > I am performing a free energy calculation using FEP by
> NAMD2.7b2. When I perform a FEP of a molecular in vacuum, I
> can get a reasonable result, eithor forward or backward.
> However, when I solvate the molecular in aqueous, it is
> surpringly that both direction yields the same results! that
> is:
> >
> > lambda     
>    forward       
>    backward
> > 0.10000        9.59495 
>          8.37989
> > 0.20000       21.96420 
>          17.85860
> > 0.30000       34.83400 
>          27.51460
> > 0.40000       47.86530 
>          37.10170
> > 0.50000       60.88580 
>          46.55320
> > 0.60000       70.98240 
>          52.97450
> > 0.70000       80.99340 
>          59.52030
> > 0.80000       90.86700 
>          67.15570
> > 0.90000      100.75300   
>        77.51150
> > 1.00000      110.93300   
>        93.77790
> >
> > I think there must be something seriously wrong in my
> configuration. Yet I cannot find out. My conf file is:
> >
> >
> #############################################################
> > ## JOB DESCRIPTION         
>                
>            
>    ##
> >
> #############################################################
> >
> > # Alchemical transformation to compute the free erengy
> of
> > # A-CH3 <-----> A-H   in a Water
> Box
> >
> >
> >
> #############################################################
> > ## ADJUSTABLE PARAMETERS       
>                
>            ##
> >
> #############################################################
> >
> > structure         
> ./a_CH3_H_aq.psf
> > coordinates       
> ./a_CH3_H_aq.pdb
> >
> > set temperature    298
> > set outputname     a_CH3_H_aq_eq
> >
> > # Important for restart a simulation
> > firsttimestep      0
> >
> >
> >
> #############################################################
> > ## SIMULATION PARAMETERS       
>                
>            ##
> >
> #############################################################
> >
> > # Input
> > paraTypeCharmm        on
> > parameters         
> ./par_all27_prot_lipid_MAB.prm
> > temperature     
>    $temperature
> >
> >
> > # Force-Field Parameters
> > exclude         
>    scaled1-4
> > 1-4scaling          1.0
> > cutoff             
> 12.
> > switching       
>    on
> > switchdist          10.
> > pairlistdist        13.5
> > LJcorrection        yes
> >
> > # Integrator Parameters
> > timestep           
> 1.0  ;# 1fs/step
> > rigidBonds          all 
> ;# needed for 1fs steps
> > nonbondedFreq       1
> > fullElectFrequency  2 
> > stepspercycle       10
> >
> >
> > # Constant Temperature Control
> > langevin           
> on    ;# do langevin dynamics
> > langevinDamping     1 
>    ;# damping coefficient (gamma) of 5/ps
> > langevinTemp        $temperature
> > langevinHydrogen    off    ;#
> don't couple langevin bath to hydrogens
> >
> >
> > # Periodic Boundary Conditions
> > cellBasisVector1    22.   
> 0.   0.
> > cellBasisVector2 
>    0.   20.   0.
> > cellBasisVector3     0. 
>   0   22.
> > cellOrigin  24.95806884765625 -19.506383895874023
> -13.278444290161133
> >
> > wrapAll         
>    on
> >
> >
> > # PME (for full-system periodic electrostatics)
> > PME             
>    yes
> > PMEGridSizeX        25
> > PMEGridSizeY        25
> > PMEGridSizeZ        25
> >
> > # Constant Pressure Control (variable volume)
> > useGroupPressure      yes ;# needed for
> rigidBonds
> > useFlexibleCell       no
> > useConstantArea       no
> >
> > langevinPiston        on
> > langevinPistonTarget  1.01325 ;#  in bar
> -> 1 atm
> > langevinPistonPeriod  100.
> > langevinPistonDecay   50.
> > langevinPistonTemp    $temperature
> >
> >
> > # Output
> > outputName         
> $outputname
> >
> > restartfreq     
>    2000     ;# 1000steps
> = every 2ps
> > dcdfreq         
>    2000     ;# 1000steps
> = every 2ps
> > xstFreq         
>    2000     ;# 1000steps
> = every 2ps
> > outputEnergies      2000 
>    ;# 1000steps = every 2ps
> > outputPressure      2000 
>    ;# 1000steps = every 2ps
> >
> >
> >
> #############################################################
> > ## EXTRA PARAMETERS         
>                
>               ##
> >
> #############################################################
> >
> >
> >
> #############################################################
> > ## EXECUTION SCRIPT         
>                
>               ##
> >
> #############################################################
> >
> >
> > alch             
>       on
> > alchType           
>     FEP
> > alchFile           
>     a_CH3_H_aq.fep
> > alchCol           
>      B
> > alchOutFile         
>    a_CH3_H_aq_forward.fepout
> > alchOutFreq         
>    10
> >
> > alchVdwLambdaEnd        1.0
> > alchElecLambdaStart     0.5
> > alchVdWShiftCoeff       6.0
> > alchDecouple           
> off
> >
> >
> > alchEquilSteps          100
> > set numSteps           
> 500
> > set windowwith          0.1
> >
> > for {set i 0.0} { $i < 1.0 } {set i [expr
> $i+$windowwith]} {
> >
> > alchLambda $i
> > alchLambda2 [expr $i+$windowwith]
> >
> > minimize 100
> > run $numSteps
> >
> > }
> >
> > Can anyone give me some suggestions? I can provide
> anything detailed. Thank you!
> >
> >
> > 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