Hello all, I am experiencing great difficulty trying to simulate a crystal-amorphous phase transition in ice. After considerable trial and error, I am hoping that someone in this forum may have some insight as to what I should do.
So the premise of the simulations is to gradually pressurize ice 1h and have it transition into HDA ice. I am doing this via the NPT ensemble, where I am keeping the temperature at 100K, and increasing the pressure from 0kbar to 14kbar. That is, I run the simulation at one pressure, allow the system to equilibrate, and then using the final state of that simulation, I start a new simulation at a greater pressure. An example fix npt command I am using for my simulations is:
fix mynpt all npt temp 100.0 100.0 $(100.0*dt) tri 9000 9000 $(1000*dt)
I suspect that this may be the problem, but I am not entirely sure. The expected result is that my ice 1h cell should collapse at around 10kbar of pressure at 100K. With my current input, I can only achieve a cell collapse when the pressure is over 30kbar. I am currently using the ātriā pressure mode as my supervisor suggested that the cell should shear a bit and the cell angles (alpha,beta,gamma) should not be fixed at 90 degrees. I have also tested the āisoā and āanisoā, but have reached the same result where the cell does not collapse when it is supposed to. Being that I am currently testing the ātriā pressure mode, I have set the tilt factors xy/xz/yz to 0 at my initial state at 0kbar, as this has been the only setting I have used that allows the cell to shear as I gradually increase the pressure. I have also verified with my supervisor that the ice 1h structure we generated for the 0kbar initial state is correct, so we know that this is not the source of this problem.
Some more important details about these simulationsā¦ I am not using the TIP4P or SPE/C water potentials, I am instead using a machine learning potential. More information of this potential can be found [here]. The LAMMPS version I am using is the August 2023 release, and I am also using the DeepMD-kit as a plugin built in with LAMMPS to interpret the machine learning potential. I cloned the most recent version of the DeepMD-kit package when compiling everything together.
If there is anything else about my input script I should mention, please let me know and I will reply with the relevant information. Any insight on how to possibly correct this issue is welcome. Thanks!
Please note that this is not so much a problem of using LAMMPS but a problem of understanding the physics of the process you ware investigating. This makes it mostly (like 99%) off-topic for a forum that focuses on LAMMPS itself.
Phase transitions in general are activated processes. They need to have sufficiently large (local) fluctuations to activate/nucleate the transition. Seeing this happen spontaneously can be rare and that makes it difficult to simulate in the straightforward way that you are attempting. The classic example is supercooled water. You can cool down water way below the freezing point and remains liquid if you handle it careful and avoid any vibrations or other handling that would result in internal density fluctuations. You can then tap on the container and the water will spontaneously freeze. This is possible even on a macroscopic level. With the length and time scales available to simulations, these effects are very much pronounced. One way to avoid this would be to do coexistance studies, i.e. where you investigate how the phase boundary changes under different conditions in a system that has already both phases. This removes the hysteresis effects.
Thus the proper way to proceed is to first learn from textbooks and suitable publications how to properly simulate, investigate, and analyze phase transitions from MD simulations and then check how to translate this knowledge to your specific case.
While having no intention to contradict Axelās statement about it not being a LAMMPS question, I just wanted to make a random comment that I found your post intelligently written ^^
Thank you for your reply. I will admit, I am just beginning my graduate studies, so there is still plenty for me to learn. I understand and appreciate the idea that this is not the place for people to tell me exactly what to do. I will, however, ask if you have any recommendations on MD textbooks in general. I will look for relevant publications about simulating phase transitions in MD as well as ice specific MD simulations in the mean time
After discussing with my supervisor again, he suggests that the transition from crystalline ice to HDA ice should be sudden, and the only agitation I should need to do is gradually increase the pressure. I will, of course, read and see if I am implementing everything properly. However, I am also testing to see if heating the system slightly will allow the system to transition at the right pressure.
i.e. where you investigate how the phase boundary changes under different conditions in a system that has already both phases. This removes the hysteresis effects.
Regarding this comment, investigating the phase boundary between HDA/LDA ice and crystalline phases will be part of this project in the later stages
You can have a look at this page for book recommendation (my personal favorite is the one by Frenkel & Smit), and this page for beginners tutorials. For ice-related content, you can also have a look at Eric Hahn personal page.
I am sure the phase transition is macroscopically sudden, but that is still no guarantee that you can observe it in a nanoscale simulation. As a heuristic: suppose you know the average lifetime of the unfavored phase when 1mm^3 of the material has been held at environmental conditions favoring the new phase. Now, 1mm^3 of material contains 10^18 cubic nanometer boxes, so in general the lifetime of 1nm^3 of the unfavored phase could be up to 10^18 times longer*. Those ā¦ are not favourable odds.
*assuming that the macroscopic phase transition immediately and irreversibly occurs once any cubic nanometer of the material has crossed the activation barrier into the favoured arrangement, so that the unfavoured phaseās lifetime is an exponential random variable, with rate constant proportional to volume.
The main problem is that I should be able to. This kind of simulation has been successful using other water potentials. For example, see figure 2 here (Mechanical instability in ice Ih. A mechanism for pressureāinduced amorphization | The Journal of Chemical Physics | AIP Publishing) where a TIP4P potential is used under very similar conditions to what I am trying to simulate. With these results, as well as the results from this simulation (https://doi.org/10.1139/cjc-2017-0201), I know that I will start with ice 1h, it will begin to form some intermediate phase as the system starts to distort (this is just a distorted version of the hexagonal system, so it isnāt a drastic change). And then at around 10.2kbar, the system will undergo a very sudden and drastic collapse to form HDA ice.
Up until this point, I have only been able to form HDA ice in my simulations by overpressurizing and then relaxing the system back to 10kbar. This kind of issue can only really mean one of two things. Either there is a problem with the potential, which I doubt would be the case. Or there is an issue with how I am using LAMMPS or building the simulations. Which is likely the case and emphasizes akohlmeyās point that I should take more time to learn about what I am doing. I am still very open to tips on how to handle ML potentials in LAMMPS simulations while I am learning
None of the papers you have linked to say that you will be able to observe the quasistatic phase transition pressure directly with your simulation methods.
Both of the papers in your most recent post describe simulations of very small systems (128 molecules for the 1992 paper; a DFT study in the Canadian 2017 paper, so likely similarly very small). Those systems are likely small enough that it is easier to overcome the nucleation barrier the smaller the system gets.
We could wade through mathematics (start with DOI 10.1063/1.1667938, p 417) but hereās a simple intuitive argument: the āhardestā part of expanding one phase into another is expanding the interphase barrier, and in smaller and smaller periodic systems the interphase barrier has less and less distance to go before it āwraps aroundā, sharply reducing the additional interphase barrier area needed to increase the volume of the favoured phase.
The more recent Deep Potential water paper does not include any of its LAMMPS scripts in either the Supplemental Material or on the GitHub repositories, but a cursory read through the SM does not show (to me) that they attempted any direct simulations of a phase transition. Rather, they calculated the free energy as a function of P and T from equilibrium simulations using Hamiltonian thermodynamic integration, and then inferred phase coexistence lines from discontinuities in the inter- and extrapolated free energy function.
In any case, as your own post indicates, you have successfully completed a hysteresis loop around the ordered-disordered transition. Such hysteresis is quite usual for statistical mechanical models of phase transition, and shows that your potential can (qualitatively) reach either phase, with the hysteresis being attributable to the transition being observed along a very rapid and thus irreversible MD simulation. This is probably as good as it gets using direct molecular dynamics methods ā you will have to look into rare event sampling to get any further.
After a lengthy amount of testing (and reading relevant material), I am still unsure if hysteresis is truly the problem. In classical potentials with smaller cells (<5k atoms), I have found no hysteresis at all, I havenāt tried running my simulations with my much larger 18k atom cell to see if there is hysteresis or not.
I have, however, verified that the potential does actually work, I just needed to fill in a knowledge gap and troubleshoot some more. Following this paper here (https://doi.org/10.3389/fphy.2017.00034), I understand now that non-classical water potentials behave differently than classical potentials. For this reason, I am unsure if hysteresis is the right term to use. Regardless, I have shown that for the larger cell, the amorphization occurs if I suddenly pressurize to 30kbar and then release the pressure. For my small cell tests, I was able to achieve amorphization at around the same pressure using gradual pressurization. Furthermore, the melting point with this particular potential is 350-370K (I am running more simulations to get a better figure for this value).
As was suggested earlier, my main issue was a knowledge gap, so I thank everyone for their feedback and suggestions of what to read into. Everything, for now, seems to be working. Just in a much different way than I had initially expected
Contacting Eric directly may be the best bet. I see that on his Youtube, there are the LAMMPS readable structure files for the ice related videos. But the custom scripts for the ice melting are not linked.
If you canāt find anything there, fortunately there will be many available examples that involve ice. I believe the DeepMD-kit code has a water/ice example in it, but that is a little outside the scope of the LAMMPS forum.
In general though, it is pretty feasible to start from scratch and play around. There is a really nice tool called Genice that you can build ice structures to use as your input structure. The only issue here is that you need to edit the output file a bit to make it LAMMPS readable. The other main thing to think about is the water potential you want to use, and this really all depends on the phenomenon you want to see. The TIP potentials are pretty good in the regard that they are pretty efficient. The main issue is that you lose out on some water interactions, and if you are looking at phenomenon that occurs outside of the thermodynamic range the model was fit in, you may see erroneous results. In either case, having a good idea of how the model works is beneficial if you want to avoid the issue I had when I first made this thread.
From there, it is mainly just defining the simulation itself. Like do you want to use NVT or NPT, what is the pressure/temperature, what calculations do you want to perform, etc.
I have contacted Eric but have not heard back yet. I am trying to use new methods to construct ice structures and be careful with the choice of potential function. Thanks for your advice.
