Re: Scaling solvent-solute intermolecular interactions

From: Jason Swails (jason.swails_at_gmail.com)
Date: Tue May 13 2014 - 16:57:14 CDT

On Tue, May 13, 2014 at 4:09 PM, Michael Bellucci <bellucci_at_mit.edu> wrote:

>
> I see...This is a good point, but then this raises another question.
> Perhaps, I am confused, and maybe you or someone else could provide
> clarity. Say for example, one wishes to compute the free energy of
> solvation of a solute A in a solvent B with thermodynamic integration. The
> mixing Hamiltonian can be written...
>
> H_lambda=(1-lambda)H_1+lambda H_2
>
> where lambda is a scale factor. Let
>
> H_1=K_A+K_B+V_AA+V_BB+V_AB
> H_2=K_A+K_B+V_AA+V_BB
>
> Then the mixing Hamiltonian can be written as...
>
> H_lambda=K_A+K_B+V_AA+V_BB+(1-lambda)V_AB
>
> which is just a standard Hamiltonian with a simple scale factor
> (1-lambda) on the solute-solvent intermolecular interaction. If the
> long-range coulomb electrostatics of solute-solvent interactions cannot be
> scaled since it is a many-bodied term as you suggest, then how exactly does
> the thermodynamic integration code in NAMD handle the solute-solvent scale
> factor? This kind of scaling is similar to what I would like to achieve.
> Thank you in advance!
>

‚ÄčIt can be scaled, just not in a naively efficient way like you would think
it could. The naive implementation of nonbonded interactions is a
sum-over-pairs approach that looks something like this in pseudo-C-code

for (int index1 = 0; index1 < numParticles-1; index1++)
    for (int index2 = index1+1; index2 < numParticles; index2++) {
        r_ij = sqrt((position[index1] - position[index2])**2); // slightly
over-simplified
        energy += charge[index1] * charge[index2] * factor / r_ij; //
factor gets the units to work out
    }
}

Clearly with this approach you _must_ have some kind of spherical cutoff if
you have an infinitely periodic system since each particle has an infinite
number of periodic images. You can take advantage of this periodicity,
though, and use a "trick" to compute the full electrostatics for this
periodic system -- the Ewald sum. The trick is to recast the problem in
two parts -- one that is computed quickly in "direct space" (i.e., via the
sum-over-pairs as shown above), and the other that converges rapidly in
frequency, or "reciprocal space" using Fourier transforms (see the Ewald
derivation and explanation here:
http://micro.stanford.edu/mediawiki/images/4/46/Ewald_notes.pdf).

As Axel said, the long-range part--since it can't be computed directly as a
sum-over-pairs--cannot be decomposed (the "frequencies" you sum over are
complex functions of all coordinates). So the naive way of computing the
pair-pair interactions between just those atoms that differ in both end
states only works in the direct space sum -- it does not work in reciprocal
space.

So what you must do instead is to compute two reciprocal sums -- one for
Hamiltonian A and one for Hamiltonian B and then mix the resulting forces
and energies you get according to lambda. You can do this because at its
root, the Ewald sum is formally equivalent to the infinite sum over all
pairs of particles and you can cast the coupled Hamiltonian in this way.

I would suggest that you look at the "pairInteraction" keyword that NAMD
has. This will compute nonbonded interactions between two groups of
particles. I had actually asked a question on this list just over a year
ago (that was never answered) about how pairInteraction worked with Ewald
sums (
http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2012-2013/2536.html).
 I did a little digging into the source code back then and found that,
indeed, NAMD allocates 3 charge grids for pairInteraction calculations
presumably so that the "exact" long-range pairwise interactions between
these two groups via energy differences between the long-range
electrostatic energy of the full charge distribution less the long-range
electrostatic energies of just the first selection and the long-range
electrostatic energies of just the second selection. This might (or might
not) do what you want it to do...

Hopefully this helps,
Jason

-- 
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher

This archive was generated by hypermail 2.1.6 : Thu Dec 31 2015 - 23:20:47 CST