From: Christopher Rowley (cnrowley_at_gmail.com)
Date: Tue Nov 04 2014 - 14:43:08 CST
Hi,
I'm trying to compute the relative free energy of Na+(aq) --> K+(aq) using
FEP. I've defined a new residue for the alchemical mutation from Na+ --> K+
The lines in my PSF file are:
2437 NA 1 S2P SOD SOD 1.000000 22.0000 0
2438 NA 1 S2P POT POT 1.000000 39.1020 0
2439 CL 1 CLA CLA CLA -1.000000 35.4500 0
The lines in my FEP pdb file are
ATOM 2437 SOD S2P N 1 -9.571 -2.575 -12.631 -1.00 0.00 NA
ATOM 2438 POT S2P N 1 -9.571 -2.575 -12.631 1.00 0.00 NA
K
Because the net charge of the ion remains constant, I want to to the FEP of
only the Lennard-Jones term. I tried to avoid the electrostatic component
with alchElecLambdaStart 1.0
The pertinent parts of the input file are:
source fep.tcl
alch on
alchType FEP
alchFile bulk.fep
alchCol O
alchOutFile forward.fepout
alchOutFreq 100
alchElecLambdaStart 1.0
alchDecouple on
alchEquilSteps 100000
set numSteps 200000
runFEP 0.0 1.0 0.0625 $numSteps
The trouble is the relative energies of these two stages seems to be
averaging to roughly zero, instead of the ~16 kcal/mol its sould be
FepEnergy: 188400 3574.2698 3574.2698 -856.1131
-855.8782 0.2349 0.0527 301.1446 -0.0545
FepEnergy: 188500 3629.3847 3629.3847 -886.2057
-886.0324 0.1732 0.0529 292.5912 -0.0542
FepEnergy: 188600 3576.6946 3576.6946 -866.7726
-866.4002 0.3724 0.0537 303.3872 -0.0537
FepEnergy: 188700 3588.1746 3588.1746 -882.1014
-881.9032 0.1982 0.0538 299.9602 -0.0536
FepEnergy: 188800 3617.2403 3617.2403 -872.7395
-872.7419 -0.0024 0.0541 295.0877 -0.0533
FepEnergy: 188900 3635.4958 3635.4958 -920.0812
-920.0259 0.0553 0.0541 297.6037 -0.0532
The final change is:
FepEnergy: 200000 3609.9986 3609.9986 -885.4448
-885.4684 -0.0236 0.1624 295.1929 0.0779
#Free energy change for lambda window [ 0.9375 1 ] is 0.0778985 ; net
change until now is 0.207182
Is there something I'm doing wrong here?
Chris
This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:22:57 CST