From: Ali Khanlarkhani (alikhanlarkhani_at_yahoo.com)
Date: Wed Dec 19 2012 - 13:06:40 CST
Dear Robert
What is matter here is the order of decoupling, it is important that you
decouple electrostatic part at first and then vdw part if you do it in the reverse order,
after decoupling vdw part, you have charged particles, and
high possibility of charge overlap.
Ali
________________________________
From: "Harris, Robert C." <rocharri_at_UTMB.EDU>
To: "namd-l_at_ks.uiuc.edu" <namd-l_at_ks.uiuc.edu>
Sent: Wednesday, December 19, 2012 8:59 PM
Subject: namd-l: charging free energies with alchdecouple
I am attempting to use the alchdecouple feature in NAMD to compute the van der Waal's and electrostatic components of the solvation energy of hexane separately with free-energy perturbation, but I think that I must have made some mistake in my input file. When I compute the van der Waal's component (by setting alchElecLambdaStart=1.0 and alchVdwLambdaEnd=1.0), everything seems fine, but the electrostatic component does not seem to be correct. The electrostatic component of this energy is computed to be about -2 kcal/mol, while I was expecting values smaller than -1.0 kcal/mol to be consistent with results in the literature. The lines where I set up the FEP are as follows:
alch on
alchType FEP
alchFile hexane_solv.fep
alchCol B
alchOutFile backward_electrostatics.fepout
alchOutFreq 1000
alchElecLambdaStart 0.00001
alchVdwLambdaEnd 0.0
alchDecouple on
alchEquilSteps 5000
set lambda 1.0
while {$lambda >= 0.04} {
alchLambda $lambda
set lambda [expr $lambda - 0.05]
alchLambda2 $lambda
run 100000
}
NAMD gives the change in electrostatic energy between lambda=1 and lambda=0.95 as -0.18 kcal/mol, but when I run a regular simulation at lambda=1, compute the electrostatic interaction energy between the hexane and water, and multiply by 0.05 I get -0.006 kcal/mol. Additionally, when I try to compute this energy with TI (by simply changing alchType to TI), the electrostatic component of du/dl was -0.11, which yields a dG of -0.11*0.05=-0.0055 kcal/mol, consistent with the second number. Another strange observation is that if I use alchElecLambdaStart=0.0, as was done in the FEP tutorial, the energy is given as infinite for all lambda windows. Did I make a mistake in my choices of alchElecLambdaStart and alchVdwLambdaEnd or in some other alchemical parameter? Thanks for any help you can give me.
Robert Harris
This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:22:51 CST