Problems with simulating phase transitions of metals

Science Talk
,
1

Dear LAMMPS community,

I am conducting simulations to investigate phase transitions in metals. I have set up a simple bulk solid of aluminum, initialized at 400 K, then heated to 1325 K and cooled back to 400 K using a temperature ramp. However, I am encountering some problems related to the setup and the simulation. The script used is as follows:

clear
units metal
dimension 3
atom_style atomic
atom_modify map array

variable tempstart equal 1325
variable tempstop equal 425
variable myseed equal 12345
variable atomrate equal 1000
variable time_step equal 0.005
variable time_eq equal 100000
variable Vol equal vol
variable Temp equal temp
variable Press equal press
variable TotEng equal etotal

variable tdamp equal "v_time_step*1000"
variable pdamp equal "v_time_step*1000"

#------------

lattice fcc 4.048914 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1

region box block 0 4 0 4 0 4 units lattice
create_box 1 box
create_atoms 1 box

pair_style meam
pair_coeff * * library.meam Al Al.meam Al

#------------

reset_timestep 0

timestep 0.003

velocity all create 400 ${myseed} mom yes rot no dist gaussian

fix equilibration all npt temp 400 400 ${tdamp} aniso 1 1 ${pdamp} drag 0.2
fix 3 all print 1000 "${Temp} ${Vol} ${Press} ${TotEng}" append vol-temp-up-down2-1000-aniso-1000-4ec-eq.txt title ""


#fix 3 all print 1000 "${Temp} ${Vol} ${Press}" append vol-temp.txt title ""

thermo 1000
thermo_style custom step pxx pyy pzz lx ly lz temp etotal vol

dump dp0 all atom 100 dump-up-down2-1000-aniso-1000-4ec-eq.lammpstrj

run 50000

undump dp0
unfix equilibration

reset_timestep 0

fix melting all npt temp 400 1325 ${tdamp} aniso 1 1 ${pdamp} drag 0.2
fix 3 all print 1000 "${Temp} ${Vol} ${Press} ${TotEng}" append vol-temp-up-down2-1000-aniso-1000-4ec-take2.txt title ""

thermo 100
thermo_style custom step pxx pyy pzz lx ly lz temp etotal vol

#dump 1 all image fcc.Al.*.cfg mass type xs ys zs id
dump dp1 all atom 100 dump-up-down2-1000-aniso-1000-4ec-take2.lammpstrj

run 167000

unfix melting

fix cooling all npt temp 1325 400 ${tdamp} aniso 1 1 ${pdamp} drag 0.2
#fix 3 all print 1000 "${Temp} ${Vol} ${Press} ${TotEng}" append vol-temp-up-down2-1000-aniso-1000-4ec.txt title ""

thermo 100
thermo_style custom step pxx pyy pzz lx ly lz temp etotal vol

#dump 1 all image fcc.Al.*.cfg mass type xs ys zs id
#dump dp2 all atom 100 dump-up-down2-1000-aniso-1000-4ec.lammpstrj

run 167000

unfix cooling

1.) The output seems to be qualitatively acceptable, but it shows some artifacts. The expected result is a hysteresis in the temperature–energy graph, due to overheating and overcooling. However, the output I obtain is the following:

Picture1
Picture11756×1053 84.3 KB

The expected behavior during a phase transition is that (e.g. during melting), the energy increases while the temperature remains constant. Here, the temperature drops rapidly while the energy mildly increases, followed by a quick increase in temperature and energy back to the temperature at which the melting started.

Picture2
Picture21756×1053 66.2 KB

2.) If I increase the number of cells to 12^3, the undesired oscillations improve:

clear
units metal
dimension 3
atom_style atomic
atom_modify map array

variable tempstart equal 1325
variable tempstop equal 425
variable myseed equal 12345
variable atomrate equal 1000
variable time_step equal 0.005
variable time_eq equal 100000
variable Vol equal vol
variable Temp equal temp
variable Press equal press
variable TotEng equal etotal

variable tdamp equal "v_time_step*1000"
variable pdamp equal "v_time_step*1000"

#------------

lattice fcc 4.048914 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1

region box block 0 12 0 12 0 12 units lattice
create_box 1 box
create_atoms 1 box

pair_style meam
pair_coeff * * library.meam Al Al.meam Al

#------------

reset_timestep 0

timestep 0.003

velocity all create 400 ${myseed} mom yes rot no dist gaussian

fix equilibration all npt temp 400 400 ${tdamp} aniso 1 1 ${pdamp} drag 0.2
fix 3 all print 1000 "${Temp} ${Vol} ${Press} ${TotEng}" append vol-temp-up-down2-1000-aniso-1000-eq.txt title ""


#fix 3 all print 1000 "${Temp} ${Vol} ${Press}" append vol-temp.txt title ""

thermo 1000
thermo_style custom step pxx pyy pzz lx ly lz temp etotal vol

dump dp0 all atom 100 dump-up-down2-1000-aniso-1000-eq.lammpstrj

run 50000

undump dp0
unfix equilibration

reset_timestep 0

fix melting all npt temp 400 1325 ${tdamp} aniso 1 1 ${pdamp} drag 0.2
fix 3 all print 1000 "${Temp} ${Vol} ${Press} ${TotEng}" append vol-temp-up-down2-1000-aniso-1000-take2.txt title ""

thermo 100
thermo_style custom step pxx pyy pzz lx ly lz temp etotal vol

#dump 1 all image fcc.Al.*.cfg mass type xs ys zs id
dump dp1 all atom 100 dump-up-down2-1000-aniso-1000-take2.lammpstrj

run 167000

unfix melting

fix cooling all npt temp 1325 400 ${tdamp} aniso 1 1 ${pdamp} drag 0.2
#fix 3 all print 1000 "${Temp} ${Vol} ${Press} ${TotEng}" append vol-temp-up-down2-1000-aniso-1000.txt title ""

thermo 100
thermo_style custom step pxx pyy pzz lx ly lz temp etotal vol

#dump 1 all image fcc.Al.*.cfg mass type xs ys zs id
dump dp2 all atom 100 dump-up-down2-1000-aniso-1000-take2.lammpstrj

run 167000

unfix cooling

but the issue mentioned in point 1 remain. Worse still, the system appears amorphous /liquid already at 400 K. Also, the hysteresis observed with 4^3 cells disappears, as the system no longer crystallizes during cooling.

In summary: What causes this unphysical setting at the outset of the simulation and at the phase transition points? Could it be due to inadequate values of tdamp and pdamp (thermostat/barostat damping parameters) ?

Picture3
Picture31710×1027 55.1 KB

Thanking you in advance, I look forward to your advice.

2

@ipetranovic This is a topic asking about the science of your simulation and not really about LAMMPS itself. Thus I will recategorize it to the “Science Talk” category.

3

@ipetranovic You are arguing as if you were studying a macroscopic system. At the time and length scales accessible to MD simulations and available computing resources, things are quite a bit different. You have significant finite size effects and the heating/cooling cycles are extremely fast. Add to that, that melting/freezing are activated processes, and you have a big mess on your hands.

I am not an expert on studying such processes (what I wrote above is just common knowledge), but you should do a careful review of the published literature. Melting/freezing a metal should not be quite as bad as e.g. water or quartz.

4

I do a lot of nonequilibrium MD and I agree with @akohlmey. You are plotting a thermal equation of state E(T), but your quench rates are very fast and likely do not allow atomic processes to eatablish well-defined pressures and temperatures along your quench trajectory.

When controlling temperature or pressure with a thermo/barostat, you need to permit time for particles to equilibriate the spatial distribution of kinetic energy and atomic volume. If you do not, then the estimators LAMMPS uses for thermodynamic quantitative will not be accurate measures of the thermodynamic limit.

If you think even slower processes are important, like phase transitions, then you may need to slow your rate even more.

5

@akohlmey @toconnor Sorry for the late response. Thank you very much for your replies.

Upon further inspection of the system, there is one more unexpected behaviour I have noticed. This simulation is made of two steps: a heating step and a cooling step. They both pipe into the same trajectory file, one after the other. However, upon inspecting the trajectory file in VMD, there is a discontinuity in the simulation at the point where the heating stops and the cooling starts. During the heating step, the system over time becomes a somewhat elongated cuboid, but on the moment of transition, the system “snaps” back to a cubic box.

I suppose I have programmed the script wrong somehow, but I don’t see it. Could you please take a look?

P.S. Sorry for the wrong category.

6

I have managed to partially find the solution to most of the problems, so I will post here in hope that this will be of use to someone.

The reason the simulation slowed down over time is the thermo command. It appears that the longer the simulation runs, the longer it takes for the thermo command to update the “live-feed” graph. Commenting out all of the thermo lines removes this problem and now the simulation proceeds at roughly the same speed throughout the run,

The reason for the discontinuity between the heating and the cooling step appears to be the fact that I have used multiple dump commands which caused LAMMPS to write multiple things into the same file at the same time.