concentration constraints

[dmebane]

Hi

I have created a cluster expansion and ensemble with four sublattices. Two of the sublattices are active. I would like to fix the concentration on one of the sublattices and fix the chemical potential on the other in a SGCMC simulation. I see in the example notebook how to create a CompositionSpace object that contains such constraints, but it's unclear how to pass that object to the sampler. Trying to parse the code but it's slow-going -- any help is appreciated!

[lbluque]

Hello @dmebane 👋

If I understand correctly, I think that your situation can be resolved without relying on the CompositionSpace class. Once you have created the Ensemble, you can set the sublattices you want to "freeze" from the two active sublattices to inactive by restricting all it's sites. Then make sure that the initial occupancy you use to start your MC sets the sites of that sublattice to the concentration you need.

This will then result in MC where only sites of the single active sublattices (for which you are setting the chemical potential) will be updated. so the sites of those sublattices will never be changed, however the configuration of the originally active sublattice that you set to be inactive will still be used in computing the CE energy.

This would be something like this:

assert len(ensemble.active_sublattices) == 2
assert len(ensemble.inactive_sublattices) == 2

# say the sublattice you want to set to inactive is index 0
ensemble.active_sublattices[0].restrict_sites(ensemble.active_sublattices[0].sites)

# This should render that sublattice as inactive now
assert len(ensemble.active_sublattices) == 1
assert assert len(ensemble.inactive_sublattices) == 3

Does this sound like it addresses your use case?

[dmebane]

Hi @lbluque thanks so much for your response! I considered restricting sites as you suggest, but I was concerned that would preclude moves on that sublattice, which still has partial occupancy. However if that’s not the case then yes, this addresses our issue . . . .

[lbluque]

Ah, I misunderstood that part. So the idea is that you want to have one of the two active sublattices to have a set chemical potential and the other to have fixed composition but still allow transitions/swaps at the composition, correct?

You're right that freezing all sites won't cut it if that's the case. I think this can be addressed by setting up table flip transitions with composition constraints as you initially suggested. @qchempku2017 may be able to help better there. (Though he is traveling currently).

A question though, are there any overlapping species in the two active sublattices, ie a species present in both?

[dmebane]

Yes, exactly -- ideally one sublattice has a fixed chemical potential and the other has a fixed composition. There is no mixing between the two.* I could fix both chemical potentials and let both compositions float, and then extract results corresponding to the compositions of interest, but that seems like it could be pretty expensive, especially because I have little sense of the right chemical potential range for the fixed-composition species.

*In fact the situation is that both species can formally be said to inhabit the same sublattice, but one occupies a set of distorted sites that change the cluster orbits. To address this we created a new sublattice consisting of the distorted sites. Experiments find that there is a site-exclusion effect, in which there is no simultaneous occupancy of both the regular site and any of the associated distorted sites. We think that we can enforce this by including a term in the CE that effectively excludes such configurations, but we will have to insert it post facto since there are no such clusters appearing in our training sets.(**) Alternatively we could exclude such clusters from the MC, if that were a possibility.

**Edit: I may have been wrong about that last part. We do seem to have a term that represents that interaction.

[dmebane]

Hi -- update: I am using a workaround which is to sweep chemical potentials of the composition I want to fix, then fitting a function to the results and fixing that function at the composition I want. This seems to work at the cost of a lot of extra computation.

[dmebane]

Hi, was just wondering if there are any updates . . . much appreciated!

[lbluque]

Hi @dmebane,

Apologies for the long period of inactivity and thanks for the ping - this task had slipped my mind (both myself and @qchempku2017 started new positions)

I'll aim to draft a PR this week or next!

[dmebane]

ok thanks so much!

[kamronald]

Hi @dmebane
Just looked into this, sorry for the delay. We have a CompositeUsher already implemented, which in principle can solve your problem. However, it is currently not very user-friendly, as you have to manually create a list of MCUshers to pass into the Sampler constructor. Also, I don't think it yet works for the "table-flip" step type.
Here is example code to make a Sampler with a CompositeUsher, assuming you have an ensemble made already:

from smol.moca import Ensemble, Sampler
from smol.moca.kernel.mcusher import mcusher_factory, available_step_types

sublattices = ensemble.processor.get_sublattices()
usher1 = mcusher_factory(usher_type='flip', sublattices=[sublattices[0]])
usher2 = mcusher_factory('swap', [sublattices[1]])
usher3 = mcusher_factory('swap', [sublattices[2]])
usher_l = [usher1, usher2, usher3]
sampler = Sampler.from_ensemble(ensemble, step_type='composite', mcushers=usher_l, temperature=1000)

I have tried making one of the ushers with the 'table-flip' step type using
usher1 = mcusher_factory('table-flip', sublattices[0])
, but am running into an error which I am trying to work out. In the meantime, if the "flip" and "swap" MCUsher types work for you, then you can go ahead to use that.

[dmebane]

Hi Ronald -- thanks so much! Just in time as I have received some new CE results and need to run some MC. So to be clear, the 'flip' sublattice will run with fixed chemical potential, while the 'swap' sublattices will run with fixed concentration?

[dmebane]

Works! Thanks so much for the assistance, everyone.

[lbluque]

Glad to hear it!

Thanks @kamronald for looking into this!

Please, let us know if you have any issues or suggestions in improving the user interface @dmebane

[dmebane]

Hi @lbluque thanks so much. Indeed there is a kind of coda here -- the method works in the manner outlined by @kamronald, but we do have to fix the concentration on the sublattice that is using the 'swap' usher. We do that using the Composition subclass from pymatgen:

from pymatgen.core.composition import Composition
sublattices = ensemble.processor.get_sublattices()
compositions = [sublattice.composition for sublattice in ensemble.sublattices]
compositions[0] = Composition('Li0.291015625')
compositions[1] = Composition('La0.5703125')

which works fine until we want to add an Ewald term. For that we have to name the species with formal charges -- 'Li+', etc. -- but Composition won't accept that, returning an error something like 'you can't use + in a concentration'. Any workaround there?