RE: Scaling solvent-solute intermolecular interactions

From: Michael Bellucci (bellucci_at_mit.edu)
Date: Wed May 14 2014 - 08:50:36 CDT


Thank you Jason. Your approach makes quite a bit more sense to me. I did see the pairInteraction keyword that NAMD has and I will look into using it further. As a quick hack, that is, without modifying the source code, I wonder if it would make sense or be possible to use the thermodynamic integration portion of the code to do the scaling. That is, would it make sense to turn on thermodynamic integration, flag the solute or solvent to be perturbed in the pdb, and then run MD with a particular value of lambda? Theoretically this should have the same effect, but as Axel suggested, accuracy may be a concern. It isn't a very satisfying solution, but it might suffice for a quick and dirty calculation.

Mike


________________________________
From: Jason Swails [jason.swails_at_gmail.com]
Sent: Tuesday, May 13, 2014 5:57 PM
To: Michael Bellucci
Cc: Axel Kohlmeyer; namd-l_at_ks.uiuc.edu
Subject: Re: namd-l: Scaling solvent-solute intermolecular interactions




On Tue, May 13, 2014 at 4:09 PM, Michael Bellucci <bellucci_at_mit.edu<mailto: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 : Wed Dec 31 2014 - 23:22:24 CST