Re: charging free energies with alchdecouple

From: Jérôme Hénin (jerome.henin_at_ibpc.fr)
Date: Thu Dec 20 2012 - 03:06:05 CST

Hi Robert,

To separate the LJ and electrostatic parts of the transformation, one would typically set

alchElecLambdaStart 0.5
alchVdwLambdaEnd 0.5

which would mean that LJ (de)coupling takes place between lambda = 0 and 0.5, and electrostatic (de)coupling happens when lambda goes from 0.5 to 1. The language might be confusing because the code is written to allow atoms to be decoupled while others are being coupled - for those two groups, the meaning of lambda is reversed. The convention in the documentation and keywords is to describe the coupling transformation (hence the ...start and ...end in the keywords mean, "begin to couple", and "finish coupling", respectively).

Cheers,
Jerome

----- Original Message -----
>
> 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