That doesn't sound right, simulations initialized with stable or metastable polymorphs are not expect to immediately undergo transitions like that. At that temperature I'd expect at best some minor jitter an no loss of the crystalline structure.

Potential issues:

1) your structure is not right and has issues either in the lattice or the hydrogen directions

2) force field parameter issues, eg typos

3) strange effects due to the simulation cell shape or size, its better to use NPT for this kind of initial equilibration, that way you know the system will relax to a correct density and shape

4) simulation parameters, eg periodic boundary conditions missing etc. Check also that SHAKE or RATTLE for rigid molecules is activated etc.

Thanks you so much for the advice Martin Fitzner

I found out that by using the NPT ensemble, I was able to eliminate the problem of losing the crystal structure during equilibration.

However, there's still a lot of distortion of the atoms even though the simulation is at 100K. I have checked the values of the force-field and made sure to activate the SHAKE to constrain the bonds and angles as well and yet the results are not satisfactory.

I have attached my dump file, data file and LAMMPS script for this simulation. Any insight will be appreciated.

Thanks in advance.

The timestep is short (0.25 fs) and the temperature is low (100 K): SHAKE is unnecessary.

Martin Klvana

As you can see in the input the angle and bond energy coefficients are 0. Hence, turning SHAKE/RATTLE off would cause random behavior, certainly not water-like. It is called *rigid* water model for a reason, and you have to enforce that rigidity in one way or another, irrespective of temperature and timestep.

Abhay Vincent

Perhaps try another way of building your initial structure. I can recommend vitroid's GenIce: https://github.com/vitroid/GenIce

If you say NPT fixed the problem, then why run NVT at all? You should only use NVT if you are 100% certain about cell size and shape. A typical approach is to start with a short NPT and then switch and produce with NVT.

Martin Fitzner, okay.

Martin Fitzner Thank you for the advice. I will try and build another structure with GenIce. I wasn't aware of the use NPT first approach but now after your explanation, I do understand.

The arrangement of the hydrogen atoms is an important part of the story. The xray offer information about the oxygen atoms framework that determines the type of the crystalline group of ice. The proton arrangement in a 3D compact structure has to follow so called ice rule (one oxygen forms covalent bonds to only two hydrogens, and between two oxygen only one hydrogen exist, which form the hydrogen bond) correspond mainly two a stable water dimer (lets call it, strong HB) or a twisted one around the hydrogen bond (lets call it, weak HB). Therefore, the different decorations of the oxygen network by hydrogen atoms produce systems with different numbers of strong and hydrogen bonds, with different stabilities.

The initial geometry of the system is not related with the interaction force field (let's suppose that you have a good one and you have no mistakes in its implementation), therefore before to start the MD simulations you have to optimize the system geometry first. When you do the equilibration of the system you have to carefully produce the initial velocities (in old school, the initial velocities were zero, but the atoms were slightly displaced from their initial equilibrium position) in order to have no drift or initial tensions inside the system.

It is recommended to start the equilibration with a small systems, get it equilibrated, create a supercell of it (with the size what you are interested) and after that start the your simulation.

You may tray to use different potentials (rigid, flexible, reactive, considering the electrons ...) as some phases can be not accurately described by some of the potentials (I thing that most of them will have this problems).

Therefor, you have to look to your problem from different perspectives.

I suggest you the page https://sites.google.com/a/ncsu.edu/cjobrien/tutorials-and-guides/working-with-water-in-lammps

Regarding the time step of 0.25 fs. It is safe for initialization of the simulation to use a small time step (you can use even a smaller one for high temperatures). Do, let say 1000 steps, save the velocities together with the particle coordinates and then restart the simulation with a larger time step. Thus, you avoid to use a small time step for long simulations (not efficient), but you avoids the initial high velocities that can perturb your system by overlapping atoms and inducing phase transitions, or even generating numerical errors .

Dear Vincent,

I looked to your system and it does not look bad to me. All the water molecules are 4 coordinated, forming hexagonal rings and the hydrogen atoms are well arranged (ice rule!). In one word, your system has no phase transition. Even at 100 K (a low temperature), the rings do not remain as in the ideal case, as the water molecules vibrate around the equilibrium positions. You need longer simulations and more structural and dynamic analysis in order to judge the quality of your ice model and the accuracy of the used force fields.