Strugling to reproduce Ti-Al phase diagram

Alloy Theoretic Automated Toolkit (ATAT)
1

I am struggling to get ECIs for Ti-Al system. For now, I am considering only the configurational entropy.

I am using quantum espresso as my DFT engine. I have checked that kinetic-energy cutoff of 80Ry and KPPRA=1000, for both hcp and fcc lattices of Ti and Al, give me accuracy of 5 significant digits in the calculation of system’s energies (I hope this is enough for the full concentration range 0-100% of Ti in Al).
For fcc lattice I have generated 64 structures, but 11 of them I marked as errors (touch */error) as checkrelax was showing the amount of relaxation to be greater than 0.1 for those.
For hcp lattice I have 52 structures and 6 of them show the amount of relaxation to be greater than 0.1.
In the maps.log file I still see "New ground states predicted" and the crossvalidation score is greater than 0.025.

Should I keep waiting for the problem to resolve it selve? I was expecting that 30-50 structures should be enough as the manual says.

PS.
To use checkrelax I had to add these lines to the runstruct_qe that came with ATAT

# Abstract the final atomic coordinates to "str_relax.out" file
cat pwscf.out | getlines -jbt "Begin final coordinates" "End final coordinates" > pw.tmp
echo "1 0 0" >  str_relax.tmp
echo "0 1 0" >> str_relax.tmp
echo "0 0 1" >> str_relax.tmp
cat pw.tmp | getlines -jaf "CELL_PARAMETERS" | head -n 3 | awk '{print $1,$2,$3}' >> str_relax.tmp
cat pw.tmp | getlines -jaf "ATOMIC_POSITIONS" | awk '{print $2,$3,$4,$1}' >> str_relax.tmp
cat str_relax.tmp | cellcvrt -f -sig=9 > str_relax.out
rm pw.tmp str_relax.tmp
2

The problem in Ti-Al (and in other systems) is that a lot of structures are relaxing away from the initial lattice. Even if you eliminate the structures with clearly large relaxations, there are still some with moderately large relaxations that introduce a lot of noise in the fit.
In addition, the metastable ground state (GS) line may be very complex. So a lot of effort might wasted on getting it right although these GS would masked by other phases anyway.
(This not really a problem with ATAT, it is an issue with the cluster expansion formalism itself…)

Both of these problems can be reduced by considering a narrower concentration range, using the -c0=?? and -c1=?? options in maps.

3

Thanks!

4

Should I generate data for all concentration ranges separately?
First 0-0.25 to get Ti(hcp) phase. Then, in a separate directory 0.2-0.3 to get Ti3Al phase. And so forth? What if I don’t know what the phases are?

5

You would typically regroup phases that share a common underlying lattice. In this case, Ti and Ti3Al are both based on hcp, so they can (and should) be done together.

The issue is when you have, say, fcc and bcc phases, the fcc ground state line, when bcc should have been stable, is often very complex and difficult to reproduce when, in fact, we don’t even need to reproduced it well because it is masked by bcc. Also, in such cases, the fcc structures tends to overrelax to bcc-like structures and this complicates the cluster expansion convergence.

In a system where you didn’t know in advance what phase are present you would notice the issue by observing that the relaxations (you can see them with the checkrelax command) are very large (about >0.1) for many structures. If you notice that this happens disproportionately in some concentration range, it is a good indication that you should restrict the concentration range to exclude the "overrelaxed" region.

6

I have generated 108 structures for fcc lattice and I am trying to get the fit in 0 - 0.5 concentration range. checkrelax shows that some of the structures are overrelaxed so I marked them as error.

0.0000 0/str_relax.out
0.0000 1/str_relax.out
0.0000 27/str_relax.out
0.0000 28/str_relax.out
0.0000 499/str_relax.out
-nan 140/str_relax.out error
0.0025 52/str_relax.out
0.0029 2/str_relax.out
0.0030 552/str_relax.out
0.0084 1907/str_relax.out
0.0089 2009/str_relax.out
0.0098 983/str_relax.out
0.0116 30/str_relax.out
0.0124 438/str_relax.out
0.0134 7724/str_relax.out
0.0151 4/str_relax.out
0.0156 9270/str_relax.out
0.0157 46/str_relax.out
0.0194 74/str_relax.out
0.0240 1894/str_relax.out
0.0264 58/str_relax.out
0.0271 59/str_relax.out
0.0283 36/str_relax.out
0.0285 802/str_relax.out
0.0288 847/str_relax.out
0.0289 3/str_relax.out
0.0295 92/str_relax.out
0.0296 472/str_relax.out
0.0300 25/str_relax.out
0.0311 1919/str_relax.out
0.0315 91/str_relax.out
0.0341 357/str_relax.out
0.0353 11/str_relax.out
0.0373 31/str_relax.out
0.0391 480/str_relax.out
0.0394 502/str_relax.out
0.0396 455/str_relax.out
0.0410 132/str_relax.out
0.0412 49/str_relax.out
0.0432 98/str_relax.out
0.0436 9267/str_relax.out
0.0444 463/str_relax.out
0.0471 103/str_relax.out
0.0478 77/str_relax.out
0.0491 16/str_relax.out
0.0499 474/str_relax.out
0.0502 123/str_relax.out
0.0514 434/str_relax.out
0.0516 10/str_relax.out
0.0524 45/str_relax.out
0.0527 9252/str_relax.out
0.0533 182/str_relax.out
0.0534 41/str_relax.out
0.0536 83/str_relax.out
0.0551 73/str_relax.out
0.0558 53/str_relax.out
0.0592 113/str_relax.out
0.0596 9044/str_relax.out
0.0604 21/str_relax.out
0.0609 7732/str_relax.out
0.0622 2007/str_relax.out
0.0632 549/str_relax.out
0.0634 454/str_relax.out
0.0644 8/str_relax.out
0.0653 457/str_relax.out
0.0657 4960/str_relax.out
0.0660 548/str_relax.out
0.0663 61/str_relax.out
0.0665 24/str_relax.out
0.0712 572/str_relax.out
0.0720 76/str_relax.out
0.0735 124/str_relax.out
0.0743 545/str_relax.out
0.0750 22/str_relax.out
0.0782 1122/str_relax.out
0.0793 178/str_relax.out
0.0893 201/str_relax.out
0.0896 47/str_relax.out
0.0909 75/str_relax.out
0.0920 84/str_relax.out
0.0979 85/str_relax.out error
0.0999 299/str_relax.out error
0.1037 99/str_relax.out error
0.1075 180/str_relax.out error
0.1115 122/str_relax.out error
0.1171 93/str_relax.out error
0.1172 100/str_relax.out error
0.1178 18/str_relax.out error
0.1201 42/str_relax.out error
0.1261 546/str_relax.out error
0.1280 35/str_relax.out error
0.1300 6/str_relax.out error
0.1303 19/str_relax.out error
0.1366 116/str_relax.out error
0.1420 443/str_relax.out error
0.1666 166/str_relax.out error
0.1804 115/str_relax.out error
0.1932 14/str_relax.out error
0.2004 10502/str_relax.out error
0.2026 82/str_relax.out error
0.2371 68/str_relax.out error
0.2909 117/str_relax.out error
0.3173 169/str_relax.out error
0.3269 9677/str_relax.out error
0.4024 37/str_relax.out error
0.4322 67/str_relax.out error
0.4510 13/str_relax.out error
0.5082 101/str_relax.out error

Nevertheless, maps.log still shows

The internal database of structures extends up to 12 atoms/unit cell, see predstr.out
Among structures of known energy, true and predicted ground states agree
No other ground states of 12 atoms/unit cell or less exist.
Concentration range used for ground state checking: [0,0.5].
Crossvalidation score: 0.0457238

I am not sure how to deal with this. Please help.

7

You log file does not seem too bad of a news…
All ground states come out ok. The cv score may be a bit large, but sometimes this is not as bad as it may seem: check if perhaps the error is concentrated on a few high energy structures. (Look at the fit.out file, 4th column) Then it’s ok.
Another possibility is that, in the process of fitting the ground state energies better, maps has increased the weight assigned to some structures (fit.out 5th column) and their cross-validation errors is also multiplied by that weight, so the cv score is a pessimistic indicator in that case.

8

I’d like to follow up with this thread on the Ti-Al system and see if I’m able to reproduce the phase diagram to a certain precision. It seems there is a lot of valuable tricks to learn through this process.

I did many trial runs with maps and phb before, but none of them gave me a satisfactory plot on the fcc-hcp boundaries.

Attached is the phb plot. I generated quite a lot of structures. fcc: 212 with 39 relaxation > 0.1 marked error. hcp: 148, with 14 marked error. In turn, there are 54 clusters for fcc and 65 clusters for hcp. The concentration ranges were both 0~1. The cv scores are both smaller than 0.03.

The convex hull:
fcc:
[img:3lug0t8q]https://astro1.panet.utoledo.edu/~terencezl/static/20160104/hull-fcc.png[/img:3lug0t8q]

hcp:
[img:3lug0t8q]https://astro1.panet.utoledo.edu/~terencezl/static/20160104/hull-hcp.png[/img:3lug0t8q]

ECI plots:
fcc:
[img:3lug0t8q]https://astro1.panet.utoledo.edu/~terencezl/static/20160104/clusinfo-fcc.png[/img:3lug0t8q]

hcp:
[img:3lug0t8q]https://astro1.panet.utoledo.edu/~terencezl/static/20160104/clusinfo-hcp.png[/img:3lug0t8q]

fcc gs:

        #         c         E     E_fit    E_diff     E_abs  E_fit_abs
str_id
0       0  0.000000 -0.033950 -0.028530 -0.005420 -4.191520  -4.186100
27      1  0.250000 -0.421471 -0.435849  0.014378 -5.673309  -5.687686
9270    2  0.333333 -0.478037 -0.456075 -0.021962 -6.094629  -6.072667
3       3  0.500000 -0.441611 -0.433046 -0.008565 -6.787716  -6.779151
28      4  0.750000 -0.295046 -0.288760 -0.006286 -7.735418  -7.729132
1       5  1.000000  0.050567  0.059501 -0.008934 -8.484073  -8.475139

hcp gs:

        #      c         E     E_fit    E_diff     E_abs  E_fit_abs
str_id
0       0  0.000  0.000000  0.016936 -0.016936 -4.157570  -4.140634
156     1  0.125 -0.205531 -0.196083 -0.009448 -4.910235  -4.900787
570     2  0.200 -0.310133 -0.302722 -0.007411 -5.343117  -5.335706
251     3  0.250 -0.366822 -0.372809  0.005987 -5.618659  -5.624646
603     4  0.400 -0.417281 -0.406034 -0.011247 -6.325679  -6.314432
212     5  0.500 -0.414702 -0.414044 -0.000658 -6.760807  -6.760149
262     6  0.750 -0.303885 -0.303446 -0.000439 -7.744257  -7.743818
1       7  1.000  0.000000 -0.007354  0.007354 -8.534640  -8.541994

Here is the phb plot:
[img:3lug0t8q]https://astro1.panet.utoledo.edu/~terencezl/static/20160104/phb.png[/img:3lug0t8q]

I initially suspected it was because the reference energy was not the same for the two lattices, and have just fixed that, using the ref_energy.out from hcp as the ref_energy.in for fcc. All above represent the newer results. One can see that the c = 0 and 1 fcc phases have none-zero energy values, but those of c = 0 and 1 hcp phases are 0.

I did notice the energy difference between the lattices at c = 0.5, 0.75 are very close, about 0.01 eV. Does it have anything to do with this?

9

Terencelz’s results seem quite well-converged. The jump in the phb plot for the 3-(6-d1) equilibrium (which I assumed means gs 3 from one lattice and gs 6 from another) show a jump at 600K. This is due to a new phase is appearing between at c=0.3 (in spin units). So above 600K there are two equilibria: one (which the code is now following) between fcc phase at c=0 and the new phase at c=0.3 (and compositions to the left as T increases). The other equilibrium is between two hcp phases (one is the ordered Ti3Al phase and the other is the new phase, which I suspect is a disorder hcp phase). To make the code follow that one, restart phb at ~600K but specify the two gs as being disordered hcp on the left (gs=-1) and the the Ti3Al phase on the right (I presume this is gs=6). Make sure to specify the chemical potential that phb had found for T=600K as the starting condition.
Note that the experimental phase diagram does show such a transition!

10

Yes, although the grouping of the name is phb-36-(d1-xxx)-(d2=yyy).out. In this case directory1 is the fcc directory, and hcp in directory2.

It looks like it should wait till 1400 K to reach an eutectic point, if the shape matches the one you provided…

[img:3532io0b]LINK_DELETED/~terencezl/static/20160104/AlTi-van-de-walle.png[/img:3532io0b]

I did this, with

phb -gs1=-1 -gs2=6 -mu=0.252168 -T=600 -dT=50 -keV -dx=1e-2 -er=10 -ltep=5e-3

but only one point got generated. So I increased mu to 0.28 (phb-d6-1.out), 0.3 (phb-d6-2.out) and 0.35 (phb-d6-3.out), and got these three:

[img:3532io0b]LINK_DELETED/~terencezl/static/20160104/phb-2.png[/img:3532io0b]

Looks like 0.35 made it. But it’s not the chemical potential that phb had found for T=600K.

It looks like the starting point of T and mu has something to do with what gets generated. I encountered the same confusion in another more successful attempt. See post below!

11

Here is a more successful attempt. With the same set of structures for both lattices, but constraining the number of pairs and triplets in the cluster sets to be equal to those you used in https://doi.org/10.48550/arXiv.cond-mat/0201511. fcc: 3 + 2, hcp: 11 + 6.

ECI plots:
fcc:
[image deleted]
hcp:
[image deleted]

They look quite similar to the values you got in 2002.

And the final phase diagram:

[image deleted]

Let me explain… and post questions along the way.

The black lines (phb-36-, phb-67-) are phb runs without problems, except phb-67-T=1850-, in red rectangle 2. Some portions of the phb-36- and phb-67-* are denser, because I decreased -dT to 10 and increased -dx to 1e-3 to test it out and look for phase transitions.

The colorful lines (phb-3d-) in red regtangle 3 are continuation at 2200 K trying to equilibrate fcc-3 and hcp-disordered, with different mu values. I was quite surprised and confused by the fact that all of them actually equilibrate just fine and go up, but from different starting points. It echoes with my replying post above; well, kind of, because in this case I don’t really have to calculate phb-3d, because the hcp-6 in phb-36- seems to have already gone through a phase transition, with a very small change in concentration.

Q1: The question is valid though, why does it depend on the provided mu value. The way I had in mind, was that if you provide a different mu that doesn’t drop the system in an equilibrium, phb simply stops with x1 = x2, but in the fig above, all of them keep growing.

Q2: Does it have anything to do with phi0, which relies on the precision of LTE and HTE, which tend to be inaccurate in the mid-T range. However, I can’t provide that to each phase in phb in the command line options.

The last two colorful lines (emc2-* and emc2c-) in red rectangle 1 are emc2 SG-canonical and canonical runs in the hcp lattice. emc2--6-r.out executes 16 runs through 1900 K to 2200 K with a step of 20 K, dropped in the hcp-6 (Ti3Al) region and goes from right to left and tracks the abrupt change in the long range order (lro). The detected x1 and x2 at transitioning mu values are scatter plotted, and the quality is not that good. emc2c-*-6-b.out executes 13 runs through x (spin) = 0 to 0.2 with a step of 0.05, dropped in the hcp-6 (Ti3Al) region and goes from bottom to top, tracking lro. The detected T1 and T2 at transitioning boundaries are connected. Both tracking attempts agree to an extent.

Q3: Now where is the eutectic point? I can perhaps somewhat see it, at the crossover point of the black-red-blue lines. There is even a tiny dent in the black line there. But I just can’t set mu correctly to drop the system at the equilibrium between hcp-disordered and hcp-6. Is it because the equilibrium region is too small, and requires higher precision? But I can’t even see an obvious break in the black lines in the first place, like the one from more clusters above!

12

I think I can partially answer Q2 and Q3.

It looks like phb does keep using lte for ground starting states, and hte for disordered state, and the command phb -gs1=-1 -gs2=6 runs the risk of having discrepancies in the mid-T range.

Here is an illustration trying to get the little branch between the disordered phase and the ordered hcp-6 phase. It is around mu = 0.3 ~ 0.4, at T = 2000 ~ 2200 K.

I chose T = 2050 K to track the mu range from the hcp-6 phase into the left disordered phase, and got

[image deleted]

We see from the sharp drop of lro that it does go through an phase transition.

However, also tracking from the disordered phase into the right hcp-6 phase, the phi curves of both do not intersect. There is quite a difference (dashed lines). This is indeed due to the mid-T starting integration point, using lte for hcp-6 and hte for the disordered. The correction can be done by emc2 integration from a low T = 300 K for hcp-6 and a high T = 5000 K for the disordered phase, and use the ending phi values as the respective starting points -phi0 for emc2 and then compare curves. This is seen as the solid lines.

[image deleted]

It is clear that hte really needs a high T to be accurate. I even tried 10000 K and 20000 K, and they brought the two solid curves in the 0.001 eV close range (not shown). The phi0 integrated from 20000 K even lead to the disordered line (blue) to go below the hcp-6 line (red) (not shown). At the same time, curve crossing is still impossible to detect, because they decline at the same rate.

So my guess is that this small branch is tracked through lro, rather than finding the intersecting point of the two phi curves, isn’t it?

terencelz (forwarded by avgjoe)