[gmx-developers] question about crystal symmetry in simulations

[Shirts, Michael R. (mrs5pt)]

Hi, all-

I'm interested in increasing the ability of Gromacs to handle crystal
symmetries better, especially for smaller molecules. The "best" (under
some definitions of best) thing to do would make it possible to simulate
asymmetric unit cell, so that only the minimum number of copies of the
atom are simulated. This sounds like a boatload of work to get the PME
right to handle rotational/mirror symmetries, along with a bookkeeping
nightmare because for a small enough system there could be multiple copies
of any atom within any reasonable cutoff, something Gromacs can't quite
handle right now. I believe this has been done by Mike Schnieders at Iowa
for his code, so it's certainly theoretically possible. For proteins the
problem isn't quite so bad since one only need to handle the symmetry,
generally not worry about multiple copies of any given atom.

As a stopgap, I was thinking about implementing force averaging. This
doesn't improve the speed any, as you still would need to run a system
large enough so there was only one copy of any give atom within cutoffs,
but does allow us to impose symmetry at least between unit cells, if we
really wanted a simulation of just the unit cell. But with small systems,
one doesn't need to be too fast.

I would probably start with just the full unit cell, as only translational
symmetry would need to be considered, and the forces on each atom would be
a simple average of the forces on every other symmetric atom.  This should
maintain symmetry if one starts with symmetry, though I suppose there
could be issue with numerical instability and drift - My instinct is that
it should self-correcting, but I could be wrong. This wouldn't really help
protein crystal simulations at all, at least until non-translational
averaging was added in, which does sound like a lot more work.

I would probably implement it by using index groups -- one would simply
indicate which index groups should be considered atoms that should be
averaged.  This would result in minimal dependencies on anything else in
the code.  

One could in fact imagine calculating only the short range forces between
one atom per index group and the other atoms, and then copy that force
over to other symmetric atoms, thus getting some of the time savings
without messing with the code too much, though you'd still need to do the
PME over all the atoms in the entire cell, and it would be a pain to
integrate with domain decomposition.

So question: is this functionality at all useful, or should I just stick
it in a side branch and plan on not incorporating this?  Am I missing
something already in GROMACS?  Someone already working on something
similar? Any other thoughts or cautions on before I dive in?

Best,
~~~~~~~~~~~~
Michael Shirts
Associate Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821