I’m trying to calculate the phase boundaries involving a metastable phase in cluster expansion (several meV/atom above the ground state convex hull). Since it’s slightly above the cluster expansion convex hull, the Monte Carlo structure relaxes away from the phase of interest when we use PHB code, judging from the output composition. Canonical Monte Carlo with EMC2 code doesn’t work either for the same reason. Similarly, I can obtain the free energy convex hull using semi-grand canonical Monte Carlo using EMC2, but not the free energy curve of each phase independently.
My understanding is that it involves sampling the unstable/metastable region of the configurational space. From previous publications [Phys Rev B. 96 (2017) 134204; Phys Rev B. 85 (2012) 184203], I assume we have to add some constraints to the equilibration process in Monte Carlo simulations. Is there a way to do this using EMC2 code?
Before I answer your question, I think I should point out what will not work.
Unfortunately, just imposing constraints on the MC run will not help you solve the problem. The constraints would be on the global composition of the simulation cell, so this would not prevent a phase separation into two phases whose composition are very different from what you want.
It is the cluster expansion itself that must be modified. Just imposing compositions bound with the -c0 and -c1 options won’t help either, because, even though the code only tried to enforce the right ground state equilibrium in the specified composition range, that equilibrium could involve structures outside of that range. This is not an error, this is intended to give the correct equilibrium that MC runs would give you.
I should also point out that in the files you sent separately, I’ve noticed some unusual things. First, the ground state line your 3rd figure has a completely different shape as from the 1st 2 figures. That seems incorrect. I am basing my answer below on figure 1 and 2 only. Second, is B19 really a superstucture of hcp here?
Now, what is the solution? You need a cluster expansion where the ground states are the metastables states you want. This involves two fixes:
the CE seems to overstabilize the hcp Mg structure. You could fix that by assigning more weight to that structure. You can use mmaps (note: 2 m’s) which has some extra features. You can specify a weights.in file. You can also set the -ig option so that it does not try to enforce ground states at all, which gives you control over the weights - this is optional.
the "B19" structure is the one that prevents the GP zone from being stable (if you want the Mg-GP equilibrium). So you could artificially set the energy of that structure to force it up until the GP is stable. You could also increase its weight to make sure it happens.
Caution: you should not change the energy of the Mg structure, because your answer depends on that. You can change the energy of the B19 because you don’t want the GP-B19 equilibrium.
The first 2 figures are obtained by fitting the DFT energies after excluding the structures that relax too much (using checkrelax code). After the exclusion of several obviously off-lattice structures, B19 becomes the ground state. The 3rd figure is obtained by fitting the structures within [0, 0.1111]. Outside this composition range, it’s not reliable.
I will try to change the weights and energy of B19, so that the GP zone becomes stable, then calculate the phase boundaries.
Just a suggestion: maybe constrain the composition to include Mg and B19. That will prevent weird ground states from appearing.
BTW the -c0,-c1 options of maps are replaced by an input file crange.in in mmaps (because it’s more general).