FEP yields the same results for forward and backward computation!

From: Jun Zhang (coolrainbow_at_yahoo.cn)
Date: Wed May 12 2010 - 22:38:44 CDT

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