# How to simulate a phase transition in Molecular Dynamics?

I am trying to simulate the phase transition for gold and silicon separately using LAMMPS. I got the melting point for gold right using the below code.

units   metal
atom_style      atomic
boundary        p p p
variable        a equal 4.0782

lattice fcc 4.0782
region  box block 0 10 0 10 0 10
create_box      1 box
create_atoms    1 box
mass    1 196.97

pair_style      eam
pair_coeff      * * Au_u3.eam

minimize        1.0e-8 1.0e-8 1000 100000
min_style       cg

timestep 0.001
velocity all create 300.0 454883 mom yes rot yes dist gaussian

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

fix 1 all press/berendsen iso 0.0 0.0 100.0
fix 2 all nvt temp 300.00 2400.00 1.0

run     10000000


However, when I apply this code to silicon, I don't get the right results. I feel like I am not understanding the code and the physics behind getting a phase transition graph using molecular dynamics. So I guess my request here is if you have recommendations of some papers or books I should read to be able to perform the simulation (phase transition in molecular dynamics). I have browsed through the book " Understanding Molecular Dynamics by Dan Frenkel and Berend Smit", but still I feel like I am missing something.

edited : the output is shown below for both gold and silicon.

• Tried different boundary conditions?
– Emil
Feb 7, 2018 at 6:17
• Show us what you get with gold and what you get with silicon. Give us more details on why do you think that your simulation has failed. Feb 7, 2018 at 9:03
• (An aside: your nvt time constant is about $10\times$ too large) Feb 7, 2018 at 9:26
• @valerio92 the output is shown above, where the in gold the phase transition is shown to be at around 1300 k which is correct and I did calculate the latent heat and was correct too. however, for silicon the melting point was not correct where it was around 2000 k and the correct value is around 1600 k.
– phyM
Feb 7, 2018 at 13:27

Metals, such as gold, undergo phase transitions with very limited nucleation barrier, whereas silicon does not.

For instance, if you have periodic boundary conditions (no surfaces) and no other defects, then silicon can be superheated very easily (at least on the typical time-scale of an MD simulation).

If your aim is to measure the melting point then the correct way is to start off by creating two phases - crystalline and amorphous - in contact with each other, and then find at what temperature the two are in equilibrium. (At lower temperatures the crystal phase will grow, at higher temperatures the amorphous phase will grow).

If, however, your aim is to simulate the nucleation event then this is a very difficult problem and a very active area of research. You could simply go to a higher temperature or fix your cell size (remove the press/berendsen fix) - that will cause a phase transition, but it won't necessarily be a realistic transition. There are much more complicated procedures using rare-event sampling methods like metadynamics [1] or seeding methods [2] for this.

• I am trying to measure the melting point. so I guess I should create two phases but can you elaborate more on how to do that?
– phyM
Feb 7, 2018 at 13:35
• @maryam It should be periodic with half of the crystal melted; NPT ensemble; run the simulation at multiple temperatures (e.g. 1000 K, 1100 K, ...). For each one, record whether the crystal grew or the amorphous phase grew. You will find a threshold temperature that separates these two classes - that's the melting point. If you need a precise measure of the melting temperature you can use an order parameter (either potential energy or a geometrical one like the centrosymmetry parameter) to measure the growth of each phase. ... If you ask a specific question I can give a specific answer. Feb 7, 2018 at 16:20
• thank you! I did what you said and the results I got is: increasing the temperature each run does not affect the system. The amorphous phase does not increase even though I tried running it at 2500 K. so I think there is something wrong with my script. I would really appreciate it if I can share the script with you and help me see where I went wrong.
– phyM
Feb 26, 2018 at 9:46