# How to get the phase transition diagram for a coexistence phase simulation?

I am running a lammps program to get the phase transition diagram. one for gold and the other for silicon. I am measuring the melting point by investigating the growth of each phase at different temperatures. by this method, I get a rough value of the melting temperature. but I want a precise measure of the melting point. I know I have to "use an order parameter to measure the growth of each phase". So my question is how to do that? should I calculate the potential at the end of each simulation for the different temperatures and plot potential energy vs temperature diagram. or is there a more accurate/effective way to do it? below is the code I used in my simulation for gold"

units   metal
boundary    p p p
atom_style  atomic

variable    a equal 4.0782
lattice fcc \$a

region box block 0 10 0 10 0 10
create_box 2 box
create_atoms 1 box

region       box2 block 0 5 0 10 0 10
set      region box2 type 2

group        solid type 1
group        melt type 2

mass * 196.97

pair_style  eam
pair_coeff  * * Au_u3.eam

timestep 0.002

neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes

thermo  100
thermo_style    custom step pe ke etotal temp vol press density atoms

dump 1 all custom 500 dump.all id type x y z
dump_modify 1 first yes

# melt just the 'melt' region first

velocity melt create 1000.0 556578
fix 1 melt nvt temp 1000.0 5000.0 0.2
run 10000
unfix 1

# now relax the entire structure to its local minimum

min_style cg
minimize 1e-4 1e-6 100 1000
run 0

# now you can test which phase melts or grows at a given temperature, e.g. 1500 K

velocity all create 1500.0 1234
fix 1 all nvt temp 1500.0 1500.0 0.2
fix 2 all press/berendsen iso 0.0 0.0 1000.0
run 5000