Si-Ge phase diagram converging point > 450 K

Hello,

I am trying to reproduce the phase diagram of Si-Ge, but numerous attempts result in the phase converging point larger than 450 K. I tested the whole workflow on two different supercomputer clusters, both the stable and 3.07 version of atat, and different tags (SIGMA 0.1, 0.05, KPPRA 1000, 2000, KPAR/NPAR) and made sure that VASP is printing out valid results, but all of them gave me a converging temperature higher than in https://arxiv.org/abs/cond-mat/0201473, which is around 300 K.

Here are my inputs:

vasp.wrap

[INCAR]
PREC = high
ISPIN = 1
ISMEAR = 0
SIGMA = 0.1
EDIFF = 1E-5
NSW = 41
IBRION = 2
ISIF = 3
EDIFFG = -0.01
KPAR = 3
NPAR = 1
KPPRA = 1000
DOSTATIC
USEPOT = PAWPBE

lat.in

6.062178 0.000000 0.000000
0.000000 6.062178 0.000000
0.000000 0.000000 6.062178
0.000000 0.500000 0.500000
0.500000 0.000000 0.500000
0.500000 0.500000 0.000000
0.000000 0.000000 0.000000 Si,Ge
0.250000 0.250000 0.250000 Si,Ge

Output:

cluster.out

1
0.000000
0

2
0.000000
1
1.000000 1.000000 1.000000

4
2.625000
2
1.000000 1.000000 1.000000
0.750000 0.750000 1.250000

12
4.286607
2
0.250000 0.250000 0.250000
-0.250000 -0.250000 0.250000

12
5.026492
2
0.250000 0.250000 0.250000
-0.500000 0.000000 0.500000

6
6.062178
2
0.250000 0.250000 0.250000
-0.750000 0.250000 0.250000

12
6.606105
2
0.250000 0.250000 0.250000
-0.500000 -0.500000 0.000000

24
7.424621
2
0.250000 0.250000 0.250000
-0.750000 -0.250000 0.750000

12
7.875000
2
0.250000 0.250000 0.250000
-1.000000 0.000000 0.000000

4
7.875000
2
0.250000 0.250000 0.250000
-0.500000 -0.500000 1.000000

12
8.573214
2
0.250000 0.250000 0.250000
-0.750000 -0.750000 0.250000

eci.out

0.043204
0.000136
-0.000657
0.000213
-0.001500
-0.001049
-0.001100
-0.000343
0.000459
-0.001566
0.000220

This is used to generate the phase diagram.

phb -gs1=0 -gs2=1 -dT=25 -k=8.617e-5 -dx=1e-3 -er=50 -ltep=5e-3 -o=ph01.out

Well, it looks like a matter of potential choice.

(PAW-)PBE gives this overestimation, compared with (PAW-)LDA, which was used by https://arxiv.org/abs/cond-mat/0201511.

Well, that’s a relief, since all curious factors are addressed.

Still regarding this case, when I was trying to run

emc2 -gs=0 -mu0=0.5 -mu1=1.5 -dmu=0.1 -T0=280 -T1=500 -dT=10 -keV -dx=1e-3 -er=50 -o=mc01.out

the calculated phi, E_lte-mu*x_lte and phi_lte are very large (3.4E38), and from the output, the calculation is painfully slow.

I tried from a small T0, 10 K to a large one, 280 K, which had resolved the problem of other systems, but this case it did not work.

The full output is

When phi_lte is equal to ~1e38 that just means that the low T approximation is used in a regime that is not low T. Solutions:

  1. start at a lower T
  2. is -gs=-1 (which will trigger the code using the high T expansion for the starting point instead)
  3. for the initial value of phi with -phi0=… (here different runs starting at different point must use the same reference, this is more tricky to achieve). This is mostlly useful to continue a stopped calculation with the last phi obtained.

Near a critical point, calculations are slow, that’s normal. Maybe use -eq=… -n=… instead of -dx=… to force finishing in a given time rather than a given precision.

a) I assume lte, hte, and mf are not cumulative on the integration path. It’s just one chemical potential and temperature point and there you have a value.

b) In this Si-Ge case, Not only lte, but Phi from the integration is also humongous.

c) I started at 10 K, but the problem still persists.

Got it. After a month of reading and trial and error (and despair). Today it’s started working.

Coupling the results of mu and x from the phb run, and the emc2 -h docs:

The starting and step value of mu should be chosen very carefully. Too large (>0), and even too large of a step (> 0.01), you easily get pass the transition point, and sometimes far enough, even to the zone where phi_lte and e_lte blow up to 1e+38.

A good operational sense of magnitude is vital. Here is the code that got me some nice free energy curves around the transition point of phase 0. -mu0 and -dmu has to be really small though. -mu1 won’t be reached, because the transition would have already taken place right above 0.

emc2 -gs=0 -mu0=-5e-3 -mu1=0.01 -dmu=5e-4 -T0=200 -T1=400 -dT=20 -keV -dx=1e-3 -er=50 -o=mc01.out

Thanks for sharing the solution! and sorry for leaving in the dark for so long (you saw the reason in the announcement heading!)