[gmx-developers] question about crystal symmetry in simulations

[Berk Hess]

Hi,

What is the purpose of such a simulation?
Putting a single copy of a small molecule in a crystal geometry will 
allow so little freedom of motion that is most cases nothing will 
happen, even when it should. Multiple copies will help here and might 
even be necessary to get a unit-cell that is large enough for a 
reasonable cut-off, although indeed LJPME might help. You probably need 
more advanced techniques, e.g MC moves, to not simply get out what you 
put in.

Additionally there are technical issues with many crystal symmetries. I 
have been asked multiple times if I could implement the rotational 
symmetry for typical multimeric membrane proteins. But this results in a 
singularity at the rotational point in the middle, which is not 
compatible with MD. It might work with energy minimization though 
(although it will require a lot of work).

Cheers,

Berk

On 10/26/2014 08:55 AM, David van der Spoel wrote:
> On 2014-10-25 23:28, Shirts, Michael R. (mrs5pt) wrote:
>> 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.
>
> What's the problem with simulating a whole unitcell? Especially for 
> small systems that won't be a problem. By using PME and LJPME with a 
> fine grid you can reduce the cutoff.
>
> One problem an old coworker of mine found was that the unit cell 
> dimensions are not maintained perfectly (but almost) which could be 
> caused by the angle/dihedral potentials not matching the geometry 
> perfectly (that means if you have a couple of angles that are 106 
> degrees, and the force field thinks it should be 108 degrees this 
> could cause high pressure and deform the unitcell somewhat).
>
> Symmetry breaking (i.e. different copies of an atom that should follow 
> crystal symmetry) will occur to some extent, but this is natural. If 
> you would implement force averaging to compensate for force field 
> problems that would not in my opinion make the simulation more 
> trustworthy. If I were to referee a paper using that I would demand 
> that the simulation was redone without averaging to test stability :).
>>
>> 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
>>
>
>