Hello all, I have recently calculated the Born matrix and am now beginning to analyze my data and see some unexpected behaviour in the matrix elements that are calculated. Here is the relevant portion of my input script for this calculation
units metal
boundary p p p
atom_style atomic
neighbor 1.0 bin
neigh_modify every 10 delay 0 check no
plugin load libdeepmd_lmp.so
read_restart restart.1895000
pair_style deepmd H2O-Phase-Diagram-model_compressed.pb
pair_coeff * *
velocity all create 100.0 23456789
compute myEp all pe/atom
compute myStress all centroid/stress/atom NULL
compute myPress all pressure NULL virial
compute myBorn all born/matrix numdiff 1.0e-5 myPress
From here, I am performing NPT simulations to monitor the phase transition from ice 1h to amorphous ice. When comparing my results to literature values, I see that the diagonal components that are calculated agree quite well (I am interested in C11, C22, C33, C44, C55, and C66). However the off-diagonal elements to not agree at all (I am interested in C12, C13, and C23 here). Additionally, when I am looking at the behaviour of C44, C55, and C13 as a function of pressure, it responds the opposite way than what is expected. Namely, I am expecting that these values decrease at this point, while they seem to increase from the data I collected.
So I am curious about a couple of things here. Firstly, is my implementation of this calculation correct? From what I have read in the documentation, and troubleshooted with a couple of errors along the way, it appears to be fine. But if someone can see something wrong please let me know. If my implementation of the calculation is correct, I am curious to know how the calculation is done. For example, is it assuming I have an orthorhombic system for some reason? Or is the calculation performed with some alternative definition?
You are not likely to get much help with running DeepMD based simulations since that is not part of LAMMPS and we don’t know much about its implementation etc. There is a large variety of machine learning potentials available in LAMMPS, perhaps one of them is worth considering as an alternative. That said, I don’t have the expertise to rank or compare those.
I suggest you set up a test system where the analytical version of computing the born matrix is available and that is similar to your water system, if possible. Then you can use the numerical differences version of the compute and see how you need to tweak the parameters to get results that are consistent with the analytical results. Then you can check if those can be transferred to other potentials, e.g. some known water potentials and still give consistent / expected results and finally redo the calculation with your DeepMD pair style and potential (or try some other ML pair style, if there are issues).
I will check out the ML potentials available and see if any will work out. I think that is a very good suggestion, I should be able to get a comparable system with an empirical model. Something like TIP4P or TIP5P should work well. Either way, I think this is a very good place to start troubleshooting. Thanks!
I see no error in your script in the invocation of the compute commands. However the results are sensitive to the statistics. Which values do you say “respond[s] the opposite way than what is expected”? Is it the Born term or the full C value including the stress fluctuations? And how do you compare with other methods? The stress fluctuation term can be long to converge for some values of the matrix depending on your potential, crystal structure and simulation setup.
The formula to get the elastic constant in the NPT ensemble are a bit different than the one in NVT and include the compressibility of the system if I’m not mistaking, see this article for more information.
So in this case, the values that are responding strangely are the C44, C55, and C13 matrix elements as a function of temperature (I have now realized I made a mistake in my first post). At the moment, what I have done is look at many of these matrix elements over time (as well as the thermo output like PVT and energy) to ensure convergence (I run for 0.25ns at least) before moving on to the next temperature point. At some point, I expect that each value should decline, but certain Born matrix elements should have a drastic decrease at a particular temperature. I am instead observing these matrix values to increase at the transition temperature
At the moment, I have not compared anything with my system. However, following the suggestion of Akohlmey, I will be testing my potential by doing the calculation with a slightly different ice system. The goal right now is to perform the exact sample calculations with an emprical water model like TIP4P or TIP5P which should have an analytical calculation available. With that, I can compare the results of with empirical with experimental results to make sure they are reasonable. And then compare the results of the ML potential I am using with those of the empirical potential to see exactly what is different and by how much.
For added information, I am currently looking at the transition between ice 1h and VHDA ice. What I will do to compare the results of the calculation is to look at the transition between ice 1h and HDA ice, as the elastic behaviour is well characterized for that transition at this point. I will also look at the paper you have attached to see how else I can approach the problem
This does not really answer my questions. The physics tells how the complete C values (Born term - Stress fluctuations) change over time. Not the individual Born values. So I wouldn’t make guesses solely based on a single one of them at finite temperature. Moreover, when mentioning the convergence, I was speaking about the convergence of the fluctuation term with the number of stress tensor terms you used, which can be quite long. I personally recommend recording the virial stress tensor at every timestep for roughly 1 million steps to compute the full correlation matrix. But this also depends on the potential and temperature as I mentioned.
Because of the long range part of TIPnP models, they don’t have analytical solutions implemented in LAMMPS as far as I know.