I have tried different ways of successfully getting a solid-liquid coexistence state for a simple Lennard Jones potential. All my approaches consist of setting up a block FCC cuboidal system, i.e. 10 x 10 x 100, and then melting half the system under anisotropic NpT while keeping the other half in the solid state. While NpT is performed, the liquid particles cannot interact with the solid particles. The liquid is then cooled to roughly the coexistence temperature based upon values obtained from literature sources. I then resample all atom velocities in both the solid and the liquid regions to match that on the approximate coexistence temperature. It should also be noted I am using a tabulated potential that described a force-shifted LJ potential where the cutoff distance for both the force and the potential energy is 2.5 sigma.
Once the NpT part of my system is run, I switch ensembles and perform my production run in NpH. The two algorithms I have tried in question are the aniso and coupled xy algorithms in NpH. As far as I know, the aniso algorithm does not work under these system conditions. The enthalpy and total energy both continually drift towards more negative values. The temperature and volume both decrease, while the system actually becomes more liquid. In long enough simulations the system, once becoming more liquid, tends to a complete solid.
The last method I have tried is NpH with coupled x and y dimensions and independent z dimension. This leads me to constant enthalpy, total energy, temperature and a constant proportion of solid to liquid. In short, this gives me a stable coexistence state.
The pressure in both cases is constant. I have also repeated simulations with the anisotropic algorithm for various sized timesteps and all give unsuccessful results.
My question is, why under these system conditions should the anisotropic algorithm fail, while coupling x and y dimensions independent of z should work? Is there something different about how the NpH algorithm is implemented in this way?
Any insight into this situation would be much appreciated.
I have done things similar to what you describe and got reasonable
results. There are lots of pitfalls, but you seem to have avoided all
of them. If you got it to work for coupled x/y, independent z, then I
would also expect aniso to work too. The only difference is in how
the 10x10 cross-sectional area fluctuates. Your description of the
main symptom is a little mysterious:
"The temperature and volume both decrease, while the system actually
becomes more liquid. In long enough simulations the system, once
becoming more liquid, tends to a complete solid."
The first sentence is consistent with a system that is initially too
hot: it melts and cools. But then in the second sentence, you say the
system cools so much that it freezes. Whatever is going on, this is
not an equilibrated system.
I suggest taking the coexisting NpH system that you generated with
coupled xy and restarting with aniso.
Also, I hope you are not doing anything strange like specifying
different values for the target stress components in different
directions. That might be a problem.
Thank you for your advice. I shall attempt some more delicate simulations along your suggested lines.
I had run simple tests though to compare the effectiveness of the anisotropic and isotropic algorithms for 256 atoms in an FCC lattice. The system was equilibrated in NpT over 400,000 iterations with a step size of 0.00025. Data from the equilibration run shows constant temperature and pressure for a the simple Lennard Jones system. The system remains as a solid. I then run NpH anisotropic/isotropic for several million iterations to see what happens. I found that for the isotropic algorithm, regardless of step size, the Enthalpy, Total Energy, Volume, Temperature and Pressure are all constant; that is they all fluctuate about a mean value. The system remains solid. However, for the anisotropic algorithms, I find that the system always tends to a liquid. The Enthalpy, Total Energy, Volume and Temperature all increase. The pressure remains constant, but the fluctuations about the mean tend to become smaller as time passes. The effects become more profound the shorter the step size.
I am not specifying different target stresses in different components, so there is no cause for worry there. The command I am using is literally “fix 1 all nph iso 1.0 1.0 0.25 mtk yes pchain 0”, where I am choosing my damping coefficient for the pressure to by 1000 times that of the step size, as recommended inside the LAMMPS documentation. Why this is happening for such simple test cases is not apparent to me. Do you, or anyone else, have any further information?