From: Floris Buelens (floris_buelens_at_yahoo.com)
Date: Sun Jun 26 2005 - 07:41:14 CDT
Hi everyone,
For the last few months I've been using NAMD to
perform alchemical FEP calculations, and have had to
deal with the infamous end-point catastrophe problem,
where at values of lambda approaching 0 or 1, linear
scaling of the energy function wreaks havoc with vdW
values when weakly interacting lambda-coupled atoms
clash with the rest of the system.
I came across the following thread (10th jan):
http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l/1694.html
This post and Jérôme Hénin's follow-up discussed the
benefits of replacing the conventional LJ function
with one that smooths out at small inter-atom
distances.
Implementing this new function for every atom-atom vdW
interaction was too much for my c++ knowledge. However
in investigating this I hit upon an alternative
approach, which I now present for comments.
I currently have a working version of NAMD where I
have added a cap that limits the vdW potential energy
value (per atom) reported to the FEP module to a given
number. This approach cuts off the extreme vdW clash
energies reported while leaving the actual energy
function untouched.
For example, in a ~40ps non-FEP simulation of my
solvated protein system the maximum vdW potential
energy value I observed for any one atom was just
below 10. In a FEP simulation at lambda ->0 or ->1 the
same values can run into the millions.
My approach now takes a configuration file parameter
vdw_cap, and any troublesome energy value exceeding
this vdw_cap (I'm arbitrarily using 20) is simply cut
off and set to vdw_cap, thereby discarding the worst
of the clashes.
I should stress that this cap is only applied to the
energy value reported, which is downstream of the
actual pairwise force and energy evaluations.
Therefore, the atom forces and the underlying motion
of the system are not affected.
I haven't tested this approach thoroughly yet but my
runs so far have shown satisfying results. Simulations
at lambda1->0 or ->1 without a vdw_cap regularly see
the vdW potential value jump up to crazy numbers for a
few dozen steps at a time, enough to force up the
dE_avg and dG reported. With my arbitrary value of
vdw_cap there were no such jumps, which will obviously
do wonders for convergence at high or low lambdas.
That's the plus side, but is this all a bit crude? Am
I losing too much information by just sweeping
shifty-looking vdW numbers under the carpet?
I look forward to comments on this approach and any
suggestions to develop the idea. If anyone's
interested in the modified source code or binaries
(linux) please email me, it'd be great to have other
people testing.
Floris Buelens
Crystallography, Birkbeck College
___________________________________________________________
Yahoo! Messenger - NEW crystal clear PC to PC calling worldwide with voicemail http://uk.messenger.yahoo.com
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:39:37 CST