Simulating molecular crystals

[aizvorski]

I'd like to simulate molecular crystals which have a fairly small unit cell (0.7nm) with PME/LJPME. What is the best way to do that?

• Have the simulation contain just one unit cell but make the PME cutoff really small to fit (eg 0.3nm) => Would any other PME parameters be needed to make this work with reasonable accuracy? Is there any theoretical minimum on PME cutoff - for example would a crazy small cutoff like 0.01nm work? (so the calculation is essentially all done in reciprocal space)
• Have the simulation contain many copies of the unit cell (eg 3x3x3) so that a normal PME cutoff (1nm) can be used, but make the copies virtual sites => Would LocalCoordinatesSite work for this, possibly with some dummy particles corresponding to the corners of the unit cell?
• Have the simulation contain many copies of the unit cell with real particles in each copy, and occasionally just replace the positions with the averaged position of all the copies to make it a perfect crystal; and/or add CustomExternalForce between the copies to keep them near-perfect all the time

The main thing I'd like to do with this type of system is minimize the whole thing, including both the positions of the molecules and the size/shape of the unit cell. Would adding a barostat (or MonteCarloAnisotropicBarostat) set to zero pressure work for adjusting the size of the unit cell?

Related #3757 and #1437

[peastman]

I don't think there's any minimum cutoff for PME. The shorter the direct space cutoff, the larger you need the reciprocal space grid to be. It might get very inefficient if the cutoff is too short, of course.

LJ interactions could be an issue. If you set the nonbonded method to LJPME, it will also use PME for the dispersion part (the $1/r^6$ term), so the cutoff won't affect its accuracy. But the repulsive $1/r^{12}$ term will continue to use a simple cutoff. If the cutoff is too short, that term will become inaccurate.

Does your application really require that every repeating unit be identical? Or would it be acceptable to simulate a larger unit cell with multiple copies of the repeating unit?

[aizvorski]

@peastman Thanks for the notes about PME/LJPME, good note about r^12. That probably mostly rules out cutoff=0.3nm for LJPME since that's about when r^12 starts becoming significant.

Is the automatic grid size okay for small cutoffs?

"Does your application really require that every repeating unit be identical?" - Essentially yes, but not required to run dynamics on identical repeating units, just to minimize them. Some number of iterations of minimize/adjust coordinates so copies are identical/reminimize would probably work.

Is there a way to minimize the box vectors as well as the positions? Or to get gradient of energy wrt box vectors? It seems that would have to be done essentially through numeric differentiation: take a small step in the box vectors (also rescaling particle positions?) and reminimize to get the new energy

[peastman]

Is the automatic grid size okay for small cutoffs?

Hopefully but I'm not certain. It uses a heuristic that I tested over a range of cutoffs, but only ones within the normal range. It might not apply well for extreme values.

Is there a way to minimize the box vectors as well as the positions?

Not directly. Calculating the gradients would require a virial, which OpenMM doesn't compute. There's also the question of defining exactly what you want to minimize. Presumably the free energy at some finite temperature?

[aizvorski]

There's also the question of defining exactly what you want to minimize. Presumably the free energy at some finite temperature?

@peastman Just potential energy to start, similar to minimizeEnergy(). Is that doable without virial?

I was hoping to get numerical gradients for energy wrt box vectors just by adding a small delta to each vector and reminimizing; just not sure how numerically stable the result would be.

Regarding virial, since that has come up for barostats before, is it hard to do that in a custom integrator?

[peastman]

The virial gives you the derivative of the energy with respect to the box vectors. It's complicated to compute analytically, involving changes to every force computation routine. For minimization purposes, though, a finite difference approximation should be fine. That requires multiple energy evaluations which would be expensive for dynamics, but not a big deal for just running a minimization.