Long Range Forces in Protein MD Simuation

Question

I am trying to understand how long range forces are calculated under the following situation.

SYSTEM-> Lets say we have a single protein surrounded by water molecules. Now I would like to run MD simulation for this system. It's NOT a crystal, just single protein and water.

Lets first consider short-range forces, i.e. Lennard-Jones. Here we first enclose (for simplicity) our single protein + water molecules in a cubic box. Lets also assume that edge of the cubic box is larger than the cutoff distance, hence the protein can not see its own image. Which is perfect.

Now lets consider long range forces. Here the primary tool discussed is Ewald Algorithm. But Ewald sum is based upon calculating potential due to images, i.e. instead of single protein, the algorithm is also taking into account the potential contributed by image proteins which is not the case in LJ potential due to cutoff. Now this algorithm is perfect for crystals, but, as per my understanding, NOT for the system I have described above.

Hence I am very confused on calculation of long range forces in my system. Can someone throw some light on this? I am using Gromacs.

The only solution I see is to calculate $\frac {1}{2}\sum_i\sum_j q_jq_j/r_{ij}$ for all pairs of the protein atoms, which will be of $O(n^2)$ complexity OR Fast Multipole Method using Tree Algorithm, and there is no other way around for my system!

Answer 1

You're almost right, but: The standard simulation methods put up with infinite images precisely because electrostatic forces are long-ranged!

As you yourself have pointed out, the standard LJ cutoffs (1-2 nm) are not acceptable for electrostatic interactions. And so, instead of using ever larger cutoffs to treat electrostatic interactions locally, we can use Ewald / grid methods to calculate these forces using FFTs of complexity $O(n \ln n)$, including long range forces from inside the same cell.

Historically MD has always been done with periodic boundaries for computational reasons: it's easier and more predictable to have a particle "wrap around" when it encounters a boundary, instead of using "open" boundary conditions which require complex algorithms for adding and removing particles at the edges and maintaining a $\mu V T$ ensemble. But the periodic boundaries dovetail nicely with Ewald / grid long range electrostatics, giving a system which is certainly mathematically consistent.

Are these models accurate? Certainly for fluids where large structures are not of concern at the molecular level. The effect of having an infinite lattice of proteins is probably small if there is enough water between them, and can be shown to be so (simply run a few runs with more water and show that the protein behaves similarly); any inaccuracy is (1) diagnosable and (2) likely far less than the discontinuities in, say, an RF treatment.

There are situations where PBC has significant effect. The best known of this is on diffusivities, where finite size corrections are known. Biophysically, when simulating lipid bilayers in 3D PBC, their effect on the membrane "ripples" must be accounted for carefully. But again, it's clear how to diagnose and treat these inaccuracies: make larger and larger systems and either hit convergence or identify asymptotic corrections. Far more cleverness is required to improve systematically on other methods for long range electrostatics.

Answer 2

The Ewald (or usually the faster particle-mesh Ewald or PME) treatment of long-range nonbonded forces seems to be standard in molecular dynamics of proteins. Typically the protein is surrounded by a "water box" (padding the protein by a certain amount, enough that it doesn't end up too close to an image of itself) and the system simulated consists of an infinite periodic grid of exact copies of this protein+waterbox. Note that typically there is no gap between copies - the water going from one copy to the next image over is continuous.

How realistic this is - that's a good question! I think the feeling of most people working in this field is that avoiding any weird effects at the edges of the system (eg very strong preference for water molecules to orient so there are no loose hydrogen bonds at the edge) is probably more important than the interaction between images. How to give a water box or water droplet an edge that behaves as though there is more water next to it, without making it periodic, seems to be an unsolved problem.

There is an alternative treatment called "reaction field" (RF) - this assumes everything outside a certain cutoff distance is a uniform polarizable dielectric (typically the cutoff needs to be a bit longer than for PME, around 2nm). It seems RF is considered less accurate than PME, and it may be (much) slower with long cutoffs as well. Note that a non-periodic water box with RF still has an edge issue - from the perspective of a water molecule at the edge, everything outside the box from distance zero to the cutoff is vacuum; only further than the cutoff it becomes the uniform dielectric.

In other words: Ewald or PME solves two issues simultaneously, a fast calculation of long range interactions and a water box with edges that behave as though they’re facing bulk water; at the cost of introducing non-physical interactions between images. The first of these could be solved in other ways, the second doesn’t really have a good alternative.

A couple of links to the documentation of various simulation packages:

• OpenMM manual - "18.6.3. Coulomb Interaction With Cutoff" describes reaction field
• Amber manual - look at "20.7.2. Particle Mesh Ewald" and "20.7.3. Using IPS for the calculation of nonbonded interactions" for an alternative - as far as I can tell, it is possible to turn off PME is Amber, but not to replace it with RF - PME-off just means vacuum.