Chemical potential in grand MC ?

[ghost]

Hi developers,

First of all, I would like to thank the developers who released the CASM code. I recently wanted to do some research on disordered rocksalt using CASM code and have the following questions.

1. I performed the grand MC simulations on a disordered rocksalt structure. I found that when I change the reference states and energies in CASM, the relation between chemical potential vs predicted stable compositions also changes. I wonder how the chemical potential is defined in CASM.

For example, if I'd like to predict the intercalation voltage of Li into LixCoO2, is it correct to calculate the voltage by using the input parameter of mu_Li in gcmc without any corrections ?

2. I found the output of the grand MC calculations from heating and cooling runs are usually different. Do you mind sharing some experience on how to interpret these difference in regards to phase stability?

3. Do we have the code for automatically tracing the phase boundaries in CASM ?

Thanks in advance.

[jic198]

@seanwong16 I am having the same problem too. My system is LixNiO2, which is very similar to your LixCoO2. I have tried using the input parameter of mu_Li in gcmc, but the voltages I got are well below the experimental values. Recalculate mu_Li by dE/dx gave me more reasonable results but I am not sure if it's the correct way to do it. I am also doing free energy integration to get the phase diagram, which involves mu_Li as well. Please let me know if you have figured this out. I am happy to discuss more.

[bpuchala]

On the chemical potentials:

The gcmc input parameter "param_chem_pot" is a parametric chemical potential param_chem_pot(a) = dg/dcomp(a), using the parametric composition comp(a), comp(b), etc. set by the composition axes selected using casm composition. The parametric composition and parametric chemical potential are used for Monte Carlo input because of the constraints caused by the fixed number of sites. To see the conversion to standard chemical potential chem_pot(A) = dg/dx(A), and the conversion between composition and parametric composition, check the casm composition -d output. For a binary alloy it will look something like:

Possible composition axes:

       KEY     ORIGIN          a     GENERAL FORMULA
       ---        ---        ---     ---
         0          B          A     A(a)B(1-a)
         1          A          B     A(1-a)B(a)

Currently selected composition axes: 0

Parametric composition:
  comp(a) = 0.5*comp_n(A)  - 0.5*(comp_n(B) - 1) 

Composition:
  comp_n(A) = 1*comp(a) 
  comp_n(B) = 1 - 1*comp(a) 

Parametric chemical potentials:
  param_chem_pot(a) = chem_pot(A) - chem_pot(B) 

The value comp_n(A) is number per unit cell of species A.

I found the output of the grand MC calculations from heating and cooling runs are usually different. Do you mind sharing some experience on how to interpret these difference in regards to phase stability?

For many systems and conditions the Monte Carlo results show hysteresis on heating / cooling or varying chemical potential based on the starting state, finite time of the calculation, temperature, etc. keeping the system from transitioning into the true equilibrium phase. The free energy equivalence determines where in the hysteresis region the actual phase boundary is.

Do we have the code for automatically tracing the phase boundaries in CASM ?

We don't release any code for this at this time.

[jic198]

@bpuchala Thanks for your reply. It seems that I cannot use the param_chem_pot for free energy integration, and I should use chem_pot instead. You mentioned chem_pot(A) = dg/dx(A), which means chem_pot is calculated from the derivative of the gibbs free energy with respect to the concentration. But from my understanding, to get the gibbs free energy, we need to know chem_pot first. So, I am at a dead-end here. I am following this paper (from Anton's group) to carry out free energy integration.

Also, according this paper (also from Anton's group), phase transition can be determined from a formation energy vs composition plot. Is it the same with the intersection of the grand canonical free energies?

Thanks for your time!

[bpuchala]

The integral given in the first paper can be carried out with composition & chemical potential as written, or parametric composition & parametric chemical potential. That follows from the definitions, the rate of change in Gibbs free energy with respect to either. From the casm composition -d output you can also make the conversion. In the example I gave,

param_chem_pot(a) = chem_pot(A) - chem_pot(B) 

I could, for example, assume that chem_pot(B) is constant, to make the conversion param_chem_pot(a) = chem_pot(A). In more complicated systems it may not be quite as simple, there could be multiplication factors also, but the principle is the same that you could pick one chemical potential to be constant, or zero for vacancies, and then solve for the others.

The phase or mixture with minimum free energy is the equilibrium phase or mixture. At 0K, this is just the formation energy. If you are unsure, please review, for example in Phase Transformations in Metals and Alloys, by Porter, Easterling, and Sherif.

[ghost]

@seanwong16 I am having the same problem too. My system is LixNiO2, which is very similar to your LixCoO2. I have tried using the input parameter of mu_Li in gcmc, but the voltages I got are well below the experimental values. Recalculate mu_Li by dE/dx gave me more reasonable results but I am not sure if it's the correct way to do it. I am also doing free energy integration to get the phase diagram, which involves mu_Li as well. Please let me know if you have figured this out. I am happy to discuss more.

Hi @jic198 , here's my experience on this problem.

In my case, the ECIs were fitted based on the calculated formation energies of LixCoO2 (x = 0 ~ 1) structures. I calculated the formation energies (E_form) from the following 2 formula:

1. E_form = E_DFT - xE_LiCoO2 - (1-x)CoO2
2. E_form = E_DFT - comp_n(Li)*E_Li - comp_n(V)*E_V - (1/2)*comp_n(O)*E_O2

I found that formula 2 gave me reasonable predicted voltage profiles, which can be derived from the equation voltage = (mu(Li) - mu(Li_0)/e). But formula 1 gave a set of ECIs with much lower CV score.

Thanks.