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

From: Felipe Merino (felmerino_at_uchile.cl)
Date: Thu May 13 2010 - 01:24:37 CDT

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 - 05:22:58 CST