[gmx-users] simulations of ice beginning to spin

[Nathan K Houtz]

Hello,

I am simulating two kinds of ice: ice Ih and ice Ic (cubic ice). Both simulations seem to have developed some rotation spontaneously and I'd like to know if I can control that somehow. I'm also simulating one in a triclinic box but the output gro file appears cubic. I'm not sure if I should be concerned about that or not. Lastly, the Ice Ic structure begins to fall apart during my simulation and I don't think it should. I'd appreciate help on any of these matters. Here are some more details:

Ice Ih is done in a triclinic box, and originally looks like this: http://imgur.com/VQg9xQz. I use a rigid TIP4P/Ice water model, constrained by shake, and simulate it in NVT for 1,000,000 time steps (2 ns) at a temperature of 217K and a density of 0.920 g/cm3. At the end, it looks like this: http://imgur.com/8uzddRH. First of all, I'm wondering if gromacs has correctly applied the triclinic box, as it appears that periodic boundary conditions have turned it into a cube. VMD is showing the box described by the .gro file. Secondly, it has been rotated, at least a few degrees, about the z axis. It did not noticeably rotate about the x or y axes. I'm not sure how to explain how the simulation would develop any angular momentum. In my .mdp file, I have gromacs get rid of angular momentum every 100 steps, but I don't think it should develop any in the first place. Any advice?

The other simulation, ice ic, didn't go quite as well. Here is the original orientation, viewed from one of the corners: http://imgur.com/i0ncpju (I realize the orthographic projections make it harder to see the actual structure, but they make it easier to see unique patterns from various angles, which I'm trying to use to determine how it has rotated). This one was also simulated for 1,000,000 timesteps (2 ns) at 217K with TIP4P/Ice constrained by shake, but at a density of 0.931 g/cm3. The resulting structure is this: http://imgur.com/Lwl8UUw. I'm not sure I've got the exact corresponding orientation in this view here, but I believe it is. It's looking down the x-axis this time. This simulation got rotated much worse. Here are the before and after views looking down the z-axis, just like I showed for the ice Ih above: http://imgur.com/a/LRQR3. With such a large number of timesteps, I cannot output enough frames in the trajectory file to view this as a "smooth" movie and identify exactly when and why these rotations happen. It's more like a very long and rapid slide show. For the ice Ic, at least, it is possible that the whole box did not rotate, but all of the molecules simply reoriented themselves. But since at least part of the box appears to still be in the ice Ic structure, I'm not sure why they would do that. Ice Ic is only metastable at any temperature, so I would assume that if it falls apart, the structure would not come back. This is another issue. I would like to still have cubic ice structure at the end of this simulation, but clearly two of the corners have disastrously broken up. However, what looks like most of the molecules appear to still retain the structure. Do you think that it could be a problem with my box size, or the orientation? I'd like to know what things I might try to prevent the disintegration of the cubic structure. 

Thanks for any help. Regards,
Nathan