Re: charging free energies with alchdecouple

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 : Mon Dec 31 2012 - 23:22:23 CST