There is something strange in calculation of the energy in a case of atom-vacancy calculations by ATAT.
maps uses number of sites rather than number of atoms for computing the energy. It looks not correct for vacancy system. I compare two systems with vacancies. Depending on the formula for the energy calculation the preference of a system changes.
The number of ground states list also looks different if I take Eformation instead of the energy from ATAT.
Normalizing by the number of sites (even if occupied by vacancies) should be work, as long as you compare energies computed in the same convention. (This convention is used in CALPHAD modeling too.)
In Monte Carlo, the chemical potential of vacancies should be set to zero.
In a nutshell: the different normalizations are just a linear transformation of each other, which cannot change the convexity of the convex hull of ground states.
First just plot two Emix for taking into account the number of atoms change and ATAT. You’ll see they are shifted along x in respect to each other.
For x>0.5 there are 2-3 more ground states (points lower than GS-line) if I plot Eform. I have no idea how to explain this issue in a publication.
And main thing. You see that A9C9 is lower by Eform and it is really true. A5C6 structure must appear at the temperatures of about 500-1000 K according to the energy difference. I do not know how to correctly compare two structures with ATAT energies. If I compare like Eform, A5C6 is lower by energy. Also phb yields me the structure A5C6 is prevailing.
And right now I must quickly decide what to do with my data. I have experimental data at me showing A9C9 structure on the scans and I have ATAT data convincing me this structure must not exist in the temperature range from 300 to 3000K.
Why the above differences appear is explained below.
In ATAT the following formula is used for the energy.
(1)
Emix-ATAT=E/NC - (1-x)E0/NC - xE1/NC
where NC is number of C atoms in the example above.
If there is a substitution to other atom, nothing really changes and the normal Emix differs only by a factor of constant.
Emix=E/N - (1-x)E0/N - xE1/N
where N is the total number of atoms in the system, in the case of substitution N equals to the number of sites.
In the case of vacancies Emix must be calculated in a different way:
(2)
Emix=E/(NA+NC*(1-x)) - (1-x)E0/(NA+NC) - xE1/(NA)
where NA and NC are the number of sites in A and C sublattices. x in the first term takes into account the change of atoms in the system.
And this formula is not more linear in respect to x. It cannot be linearly converted into the expression from ATAT.
In understand, the energy which in Eq.2 is attributed to atoms you attribute it to vacancies in Eq.1. May be there is a strong physics and deep explanation, but I have come across to a controversy in comparison of these two energies.
Are A5C6 and A9C9 different lattices (lat.in) or different structures (str.out)?
If they are different lattices, the problem may have a different origin:
The formation energies in ATAT are relative to end members on the same lattice.
If the lattice differ, you should use the file ref_energy.in to force using the same reference in different lattices.
Perhaps the issue is not in the formation energy but rather in the calculation of the composition.
You may need to calculate your own composition variable (ignoring vacancies) to be able to compare convex hulls on different lattices.
You can also compare "grand canonical energies" i.e. E-sum_i mu_i*x_i where i is an element.
mu_Vac should be set to zero.
It is the concentration from gs.out
By the way, it is strange and inconvenient that x differs for gs.out and for output of emc2/phb.
They are different lattices. I compare them with -d options in phb.
It is evident that Emix for two different systems cannot be compared directly. I use Eformation for comparison in one plot.
I cannot do this for MC calculations. Well, I can set correct value in ref_energy.in for x=1 (A5Vac6 and A9Vac9) since they have similar composition. But I have no idea which value to set there for x=0 (A5C6 and A9C9) because they have different compositions.
And this also does not solve the problem with different ground states list in the per-site and per-atom energy models.
That is what I do with y variable in the above message.
This is close to what I called Eformation above. I do not see where mu_i is set up in ATAT, do you propose to take these values somewhere? There is E-mu*x energies in emc2 output. Do you imply them? Anyway it will not help in the case of phb.
Since you have multiple sublattices in your system it’s probably easier to use the multicomponent/multisublattice versions of the commands:
maps -> mmaps
emc2 -> memc2
(sorry, there is no mphb, as it’s quite difficult to implement in full generality)
Then you can specify energies for each type of atoms in ref_energy.in in the same order as in atoms.out file (output by mmaps). Rerun mmaps.
The same ref_energy.in should work across lattices.
Then you can run memc2. There you can specify the chemical potentials in the control.in file. (see memc2 -h for help). In particular, you can set mu_vac=0 .
The composition output of memc2 should be comparable across lattices. The -g2c option lets you output free energies (instead of F-mu*x), which is more familiar to most.
You can also compare those free energies across lattices also.
Setting the energy of a vacancy to 0 in ref_energy.in is probably what you want.
Vac is the absence of C, and that has an effect on the energy that should be considered in the formation energy.
In any case, any choice of reference in the ref_energy.in file will give you the same phase diagram. It’s just a reference that shifts how you measure the chemical potential of each species.
If you are concerned that the code reports the energy per site rather than per atom, then you can easily just divide each energy by 0.5+x_C after the fact.
Just to re-emphasize an important point:
If you have a system described by the formula AB_xVac_(1-x) then
plotting E= energy of AB_xVac_(1-x) per formula unit AB_xVac_(1-x) as a function of x OR
plotting AB_xVac_(1-x) per atom (E/(1+x)) as a function of x/(1+x)
Will predict the same ground states. It may change the value of the formation energies (due to the different normalization) but it will not change whether a structure is stable relative to others.
In other words, the mapping
(E,x) -> (E/(1+x),x/(1+x))
preserves straight lines an thus curvature signs and thus stability.
I am following your suggestions to do mmaps calculations. Vacancy is included in my lat.in.
3.33 3.33 3.33 60 60 60
1 0 0
0 1 0
0 0 1
0.00000 0.00000 0.00000 A,B
0.50000 0.50000 0.50000 C,Vac
I set in ref_energy.in the formation energy of vacancy to be 0.
I found that mmaps finds the energy of formation from sites. Vac occupies a site, however it is not an atom that should be considered in the calculation of mixing energy.
If I do not specify the ref_energy.in file and run ‘mmaps’, the ‘atoms.out’ and ‘ref_energy.out’ files contain:
B -6.47738
O -11.3158
A -6.33684
Vac -1.49847
It looks very odd that the energy of Vac is -1.49847 not zero, do I need to set it to zero to get an accurate Ef? By the way, if I want to specify my own ref_energy.in file, the energy of O should be the total energy of O2 molecule division by 2 (i.e. Etot(O2)/2). And for A or B, the energy in ref_energy.in should be the total energy of A or B division by the number of the atoms in the calculation cell (Etot(A)/N(A atoms)). Is this a correct understanding to get accurate Ef of enumerated structures?