Hi there,
I am trying to determine the melting temperature of Li using the two phase coexistence method, which consists of the following two steps:
1. NVT simulation - starting from a 5x5x10 solid supercell of BCC Li (with the experimental lattice parameter of Li at RT), freezing the top half, and allowing the bottom half melt completely at ~1000 K (using selective dynamics).
2. NPH simulation - the POSCAR file is the half-solid half-liquid unit cell obtained at step (1), with the information on atomic positions only (the information about velocities is removed). Then, I set the target pressure (e.g. PSTRESS = 0) and the starting temperature (TEBEG). There is also an ICONST file included, where the unit cell is constrained to maintain the orthorhombic shape with the ratio a : b : c = 1 : 1 : 2 between its lattice parameters.
When I track the temperature during the NPH simulation, I notice that only the temperature of the first time-step is indeed equal to TEBEG, but then it suddenly jumps by 50-100K (see attached figure). This happens even if TEBEG is set to 100K, for example.
What I am looking for in this process is to test whether the temperature is monotonously increasing (solidification) or decreasing (melting), but it is hard to say so explicitly since the starting temperature is far from being the one I set.
I am aware this is a relatively small amount of atoms (500 in total, 250 in each phase), but according to different publications it seems to be possible to predict the melting temperature with such a small system using this method.
Would appreciate any insight or suggestion you might have.
Thanks!