remake dopedxpy-sc-fermi

[alexsquires]

• rename bulk_vasprun, not clear that this needs to be for a high quality dos
• use _format_defect_name() in examples (is this still a private function? If it's externally useful, should it be changed?).
• remove ..._and_save() from FermiSolver methods
• defect_levels argument when annealing becomes a boolean
• docstrings
• multiprocessings
• tests
• suppress noisy warnings
• tutorials
• grid-searching
• effective dopant response
• "exceptions" for FermiSolverPyScFermi

[kavanase]

Looking good @alexsquires. I like the mix-in/abstract class approach.

The _interpolate_chempots() function is quite nice, and I think will likely be one of the mostly-used features from these additions. From what we were discussing before in the #doped Slack, would it be possible to also integrate a grid-scanning approach too? As in, the current function linearly interpolates between two points (chem pot limits), which is nice and v useful in many cases, but if for instance we have a 3D chemical potential space, our min/max of the property of interest could also lie somewhere in the middle of this space, rather than along any linear path between two vertices/limits.

Ofc the user could iterate over all possible vertex/limit combinations if they knew what they were doing, which would be a good stab at this and get you most of the way there I imagine (kind of like a 'band structure path' approach), but still is only covering certain 'high-symmetry' paths in chem pot space. So if a grid approach was also possible (where e.g. the user can just set the chempot spacing in eV, or total number of grid points), that would also be v nice I think – should be doable with some of the scipy grid space interpolation functions I think? And/or pymatgen chemical potential diagram methods maybe?
If it was possible to implement relatively painlessly, could we add these two options? ☝️ 🙏

Sorry to be piling on the requests, but I guess the other main use case I'd imagine for this part of the code would be for the more complex defect thermodynamics analysis that py-sc-fermi allows, where you want to fix certain defect/species etc concentrations and perform some constrained solutions. Is it possible to integrate this to the fermi_solver code?

Also just fyi, about your earlier _format_defect_name() question, yes it is actually quite externally useful and I'll make it a public function in the next (minor) release.

[alexsquires]

Yeah, dw this is all on the list, should be ready to go by the end of the week, I sort of trashed everything to start again (for about the 6th time now), and just adding everything back in piecewise to make sure it's all consistent with the latest thermodynamics object.

…

On Thu, 15 Feb 2024 at 04:58, Seán Kavanagh ***@***.***> wrote: Looking good @alexsquires <https://github.com/alexsquires>. I like the mix-in/abstract class approach. The _interpolate_chempots() function is quite nice, and I think will likely be one of the mostly-used features from these additions. From what we were discussing before in the #doped Slack, would it be possible to also integrate a grid-scanning approach too? As in, the current function linearly interpolates between two points (chem pot limits), which is nice and v useful in many cases, but if for instance we have a 3D chemical potential space, our min/max of the property of interest could also lie somewhere in the middle of this space, rather than along any linear path between two vertices/limits. Ofc the user could iterate over all possible vertex/limit combinations if they knew what they were doing, which would be a good stab at this and get you most of the way there I imagine (kind of like a 'band structure path' approach), but still is only covering certain 'high-symmetry' paths in chem pot space. So if a grid approach was also possible (where e.g. the user can just set the chempot spacing in eV, or total number of grid points), that would also be v nice I think – should be doable with some of the scipy grid space interpolation functions I think? And/or pymatgen chemical potential diagram methods maybe? If it was possible to implement relatively painlessly, could we add these two options? ☝️ 🙏 Sorry to be piling on the requests, but I guess the other main use case I'd imagine for this part of the code would be for the more complex defect thermodynamics analysis that py-sc-fermi allows, where you want to fix certain defect/species etc concentrations and perform some constrained solutions. Is it possible to integrate this to the fermi_solver code? Also just fyi, about your earlier _format_defect_name() question, yes it is actually quite externally useful and I'll make it a public function in the next (minor) release. — Reply to this email directly, view it on GitHub <#46 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AFMK74XGZEUFF3TPF7MZHH3YTXZ5JAVCNFSM6AAAAABBVDPFSSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNBVHE2TIMRVHE> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[alexsquires]

@adair-nicolson , could you re-parse your Cu2SiSe3 data now that everything is relatively final and I can add the grid searching stuff back in?

[alexsquires]

@kavanase proposing this as the full functionality - do you want to take a look at the example ipynb while I finish off the tests before I get to comfortable?

[kavanase]

@alexsquires trying to blast through this and get it merged and released asap.
I'll still have a few to-dos after (like supporting per-charge outputs as in the DefectThermodynamics methods, rather than always grouped for each defect, if possible). Currently trying to beef up the tests a bit, and test actual numbers to catch cases where it's not working as expected.
One issue I've hit is with effective_dopant_concentration with the py-sc-fermi backend. I'm not sure if it's an issue with the FermiSolver(backed="py-sc-fermi") code in doped, or py-sc-fermi – would you be able to have a look at this?

/Users/kavanase/miniconda3/bin/python /Applications/PyCharm.app/Contents/plugins/python-ce/helpers/pycharm/_jb_pytest_runner.py --target test_fermisolver.py::TestFermiSolverWithLoadedData.test_pseudo_equilibrium_solve_doped_backend -- --last-failed -vv 
Testing started at 17:58 ...
Launching pytest with arguments --last-failed -vv test_fermisolver.py::TestFermiSolverWithLoadedData::test_pseudo_equilibrium_solve_doped_backend --no-header --no-summary -q in /Users/kavanase/Packages/doped/tests

============================= test session starts ==============================
collecting ... collected 1 item
run-last-failure: rerun previous 1 failure

test_fermisolver.py::TestFermiSolverWithLoadedData::test_pseudo_equilibrium_solve_doped_backend 

======================== 1 failed, 1 warning in 12.67s =========================
FAILED [100%]Testing with doped backend
Testing with py-sc-fermi backend

tests/test_fermisolver.py:552 (TestFermiSolverWithLoadedData.test_pseudo_equilibrium_solve_doped_backend)
self = <test_fermisolver.TestFermiSolverWithLoadedData testMethod=test_pseudo_equilibrium_solve_doped_backend>
backend = 'py-sc-fermi'

    @parameterize_backend()
    def test_pseudo_equilibrium_solve_doped_backend(self, backend):
        """
        Test ``pseudo_equilibrium_solve`` method for both backends.
        """
        solver = self.solver_doped if backend == "doped" else self.solver_py_sc_fermi
        single_chempot_dict, el_refs = solver._get_single_chempot_dict(limit="Te-rich")
    
        # Call the method
>       concentrations = solver.pseudo_equilibrium_solve(
            annealing_temperature=800,
            single_chempot_dict=single_chempot_dict,
            el_refs=el_refs,
            quenched_temperature=300,
            effective_dopant_concentration=1e16,
            append_chempots=True,
        )

test_fermisolver.py:562: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
../doped/thermodynamics.py:5047: in pseudo_equilibrium_solve
    conc_dict = defect_system.concentration_dict()
../../../miniconda3/lib/python3.10/site-packages/py_sc_fermi/defect_system.py:380: in concentration_dict
    e_fermi = self.get_sc_fermi()[0]
../../../miniconda3/lib/python3.10/site-packages/py_sc_fermi/defect_system.py:221: in get_sc_fermi
    q_tot = self.q_tot(e_fermi=e_fermi)
../../../miniconda3/lib/python3.10/site-packages/py_sc_fermi/defect_system.py:329: in q_tot
    lhs_def, rhs_def = self.total_defect_charge_contributions(e_fermi)
../../../miniconda3/lib/python3.10/site-packages/py_sc_fermi/defect_system.py:307: in total_defect_charge_contributions
    [
../../../miniconda3/lib/python3.10/site-packages/py_sc_fermi/defect_system.py:308: in <listcomp>
    ds.defect_charge_contributions(e_fermi, self.temperature)
../../../miniconda3/lib/python3.10/site-packages/py_sc_fermi/defect_species.py:416: in defect_charge_contributions
    for q, concd in self.charge_state_concentrations(e_fermi, temperature).items():
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = 
Dopant, nsites=1
fixed [c] = 7.004017649987961e-07
  q=+1.0, [c]=7.004017649987961e-07, deg=1

e_fermi = -0.7575499999999993, temperature = 300

    def charge_state_concentrations(
        self, e_fermi: float, temperature: float
    ) -> Dict[int, float]:
        """at a given Fermi energy and temperature, calculate the concentrations
        of the different ``DefectChargeStates`` of this ``DefectSpecies``.
    
        Args:
            e_fermi (float): Fermi energy
            temperature (float): temperature
    
        Returns:
            Dict[int, float]: key-value pairs of charge of each
            ``DefectChargeState`` and the concentration of the
            ``DefectChargeState`` with that charge, i.e.
            {``DefectChargeState.charge``: concentration}
        """
    
        var_concs = self.variable_conc_charge_states()
        fixed_concs = self.fixed_conc_charge_states()
    
        cs_concentrations = {
            cs.charge: cs.get_concentration(e_fermi, temperature) * self.nsites
            for cs in var_concs.values()
        }
        for q, cs in fixed_concs.items():
            cs_concentrations[q] = cs.get_concentration(e_fermi, temperature)
    
        if self.fixed_concentration is not None:
            fixed_conc_chg_states = sum(
                [
                    c
                    for q, c in cs_concentrations.items()
                    if q in self.fixed_conc_charge_states()
                ]
            )
            variable_conc_chg_states = sum(
                [
                    c
                    for q, c in cs_concentrations.items()
                    if q in self.variable_conc_charge_states()
                ]
            )
            constrained_conc = self.fixed_concentration - fixed_conc_chg_states
>           scaling = constrained_conc / variable_conc_chg_states
E           ZeroDivisionError: float division by zero

../../../miniconda3/lib/python3.10/site-packages/py_sc_fermi/defect_species.py:391: ZeroDivisionError

Process finished with exit code 1

[alexsquires]

It came to me in a dream.

When the “Dopant” was designated as a fixed-concentration DefectSpecies as part of the "annealing", all its associated DefectChargeState objects were also set to fixed concentrations. During the py-sc-fermi calculations, the code attempts to handle variable versus fixed charge states by determining the relative ratios of each DefectChargeState for a given defect species. However, if every charge state within that species is fixed, there are no variable states to compare against—leading to an edge case in the solver. I've fixed this in the logic for generating annealed defect systems and I will add a note to the py-sc-fermi codebase to catch scenarios where an entire defect species is fixed, so that similar cases don’t cause confusion in the future.