# Solid-Liquid coexistence under NPH

Dear Lammps users,

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.

Thanks,
Mike

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.