Stillinger Weber melting

Hello,

I am using Stillinger Weber potential to see the melting temperature of c-Si.
Stillinger Weber potential is developed in a way that it fits really well to the experimental value of melting temperature and density of c-Si.

I am using the Si.sw provided in Lammps and I am using metal units in my input file.
Basically, what I am doing is that I am ramping up the temperature from 200K to 3000K in 50,000 timesteps (dt = 0.5 fs) in NPT at 0 bar (also tried NVT). To figure out where the melting point is, I plotted T vs E_total and when there’s a “vertical” increase in E_total, that’s where I define my melting point. However, the melting temperature came out to be 2500K which is significantly higher than reported value of ~1600-1700K.

If you could give me some tips on how to improve this, I’d really appreciate that.

Thanks,

Jaeyun

Hello,

I am using Stillinger Weber potential to see the melting temperature of
c-Si.
Stillinger Weber potential is developed in a way that it fits really well to
the experimental value of melting temperature and density of c-Si.

I am using the Si.sw provided in Lammps and I am using metal units in my
input file.
Basically, what I am doing is that I am ramping up the temperature from 200K
to 3000K in 50,000 timesteps (dt = 0.5 fs) in NPT at 0 bar (also tried NVT).
To figure out where the melting point is, I plotted T vs E_total and when
there's a "vertical" increase in E_total, that's where I define my melting
point. However, the melting temperature came out to be 2500K which is
significantly higher than reported value of ~1600-1700K.

If you could give me some tips on how to improve this, I'd really appreciate
that.

please note that determining the melting point accurately in atomistic
simulations, is not so simple.
there will be a hysteresis, finite size effects, dependency of the
rate of temperature change and the type of thermostat (nose-hoover is
not very effective to move large amounts of kinetic energy quickly).

the way people typically determine the melting point is by doing a
coexistence simulations, i.e. you prepare a half molten system and
monitor whether the liquid or solid part grows at different
temperatures. the temperature where both parts remain the same size,
is the melting point.

axel.

Hi,

Thanks a lot for the help.

However, I am little bit confused. When I calculate the pair distribution function of my system at 2000K and at 2600 K, I get distinctive peaks (like a delta function impulse) for my systems at 2000K and less distinctive and broad peaks for my systems at 2600K, which is a characteristic of liquid state. Also, by using the VMD, I could also visualize that at 2600K particles were "free particles" meaning the distance each atom moves seems quite significant. These observations show that my system is in liquid state and in solid state at 2600K and 2000K, respectively. With these, I am little confused why at 2000K it acts like a solid. Any comments or tips will be very helpful.

BTW I used both NPT (allowing the system to expand at high temperature with 0K lattice constant initially) and NVT.

Thanks,

Jaeyun