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

First of all, you should choose the right order parameter. For freezing or melting transitions people often use the so-called bond order parameters: quantities that probe the structural local ordering. The Steinhardt-Nelson order parameters or related quantities are common choices.

Once you have the right order parameter, you can look at your simulation and tell how many particles are in each phase. Now you have two choices:

  • You use the order parameter to bias your simulations to directly extract the free energy of your system and, through techniques such as WHAM or its simpler sibling, histogram reweighting, which allows you to reconstruct the phase diagram.
  • You brute force the problem by performing a lot of simulations under different thermodynamic conditions (temperature, pressure, etc.). The initial conditions of each simulation are the same: you prepare a box that is split in two, with half of it filled with one phase and the other half with the other phase. You run the simulation and monitor the fraction of particles in each phase: the thermodynamic conditions for which this quantity is 0.5 determine the transition locus. This method is called "direct coexistence".

In your case I would recommend the second option, as it is the most suitable for molecular dynamics simulations. You can find a lot of references by looking for "direct coexistence melting temperature" on scientific databases. Unfortunately I am not a LAMMPS user and thus I cannot really help you with that.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.