Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Hubbard: add useful utility functions #996

Open
wants to merge 6 commits into
base: main
Choose a base branch
from

Conversation

bastonero
Copy link
Collaborator

@bastonero bastonero commented Dec 13, 2023

Useful utility functions for HubbardStructureData are added:

  • Initialize parameters via nearest neighbours analysis (using pymatgen modules for NN analysis)
  • Get a common radius which defines the nearest neighbours for all onsite Hubbard atoms

This helps issues like #991

Add the possibility of finding a shell radius defining a sphere
which contains the first (intersites) neighbours. This is helpful
when using the `hp.x` code and having a structure containing
onsite Hubbard atoms with different number of neighbours.
The Hubbard parameters can now be initialized simply
by perfoming a nearest neighbour analysis (pymatgen).
Standardization of the lattice and the atomic positions
can also be performed to have all the atoms within the
unit cell (needed for the `pw.x` routine).
Pymatgen does not order the atoms by distance.
There is not apparently any such option in the API.
Folding is then now sepatated from standardization
(the latter could produce larger structures).
@bastonero
Copy link
Collaborator Author

For the docs probably we just need to add pymatgen to conf.py

The refolding wasn't folding numbers
close numerically to 1, e.g. 0.99999999.
To circumvent the issue, a threshold is set,
and it is possible to control it via input.
Also, now the inputs for nearest neighbours
are not passed via kwargs, but via an explicit
dictionary input, which is then given via kwargs
to the NN class.
The url for pymatgen to generate interlinks
across API packages is added.
@bastonero
Copy link
Collaborator Author

@mbercx @t-reents @sphuber good to have your feedbacks here.

Example of usage:

from aiida.orm import *
from ase.io import read
from aiida_quantumespresso.utils.hubbard import initialize_hubbard_parameters

atoms = read('./input.pwi') # anything ASE can parse
s = StructureData(ase=atoms)
hs = initialize_hubbard_parameters( # returns an initialized HubbardStructureData using nearest neighbours analysis
    structure=s,
    pairs={
# onsite manifold, onsite value, intersite value, dict of intersites {kind name: manifold}
        'Mn': ['3d', 4.5, 1.0e-5, {'O':'2p', 'S':'2p'}],  
        'Fe1': ['3d', 5.0, 1.0e-5, {'O1':'2p', 'S':'2p'}],  
        ...
    },
)

Cristiano and myself used it for different calculations, and the logics seem now robust. Probably not the best input API for the initialization of Hubbard parameters, but it is the price to avoid reinventing other classes which makes even more complicated to use this method, which should be easy to use and should hide the underlying complexity of the HubbardStructureData.

"""
rmin = 0
pymat = self.hubbard_structure.get_pymatgen_structure()
_, neigh_indices, _, distances = pymat.get_neighbor_list(sites=[pymat[onsite_index]], r=radius_max)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

'get_pymatgen_structure' , could be problematic if the HubbardStructureData is 1D, or 2D

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am start wondering whether we should change the aiida-core implementation of get_pymatgen_structure instead. @sphuber, @mbercx opinions?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is wrong with StructureData.get_pymatgen_structure? Is the implementation just wrong for non-periodic structures? If so, why not just fix the bug? Or is the interface fundamentally broken and to fix it, backwards-compatibility needs to be broken? If the latter, then maybe instead of breaking it, we could look into adding a correct implementation instead. On that note, @mikibonacci has been working on a new improved implementation of StructureData that addresses its current shortcomings with respect to missing properties, such as magnetic moments etc. He is developing it here: https://github.com/aiidateam/aiida-atomistic It is still under development, but it would be nice to start thinking on adding classmethods to initialize from other types (e.g. from_ase(atoms: ase.Atoms) and from_pymatgen(structure: pymatgen.core.Structure) as well as instance methods to convert the StructureDatato other types, e.g.to_ase() -> ase.Atoms` etc.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, having a look into the get_pymatgen_structure there is a strict check. This means that any pbc different from 3D would raise error. But as far as I saw in pymatgen core modules, they actually allow for different pbc (e.g. the lattice class allows any pbc).

So one should define a Lattice instance instead of passing directly the self.cell, which is simply array like. The Lattice class also provides pbc in input, so it should be straightforward.

I don't know about older version of pymatgen, maybe they didn't have this class before? Or maybe it wasn't realised at the time?

I would simply vote to do this fix in aiida-core, since this method is used also in other codes (e.g. aiida-vibroscopy, and who know what else), and it wouldn't make sense to patch this all over the places, since then also who knows what happens if settings pbc=3D in pymatgen for a 2D system.

@AndresOrtegaGuerrero we can simply open a simple PR on aiida-core and have the fix there, and you can test it. If you understood, this should be a very easy fix, so you can try it yourself as well.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright, just opened the PR on aiida-core. @AndresOrtegaGuerrero it would be good if you tried to pull the new changes of aiida-core from this PR, and see if it solves also all the other issues raised in the other PR #998 of aiida-quantumespresso you opened.

distances = distances[sort]

count = 0
for i in range(len(neigh_indices)): # pylint: disable=consider-using-enumerate
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why not using enumerate here ?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

probably it just seemed clearer to me afterwards, but definitely possible to use enumerate

:return: (radius min +thr, radius max -thr) defining the shells containing only the first neighbours
"""
rmin = 0
pymat = self.hubbard_structure.get_pymatgen_structure()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This PR will be limited to have a new release of aiida-core where the modification you did for StructureData is merged. I would consider to do what @mbercx did in this PR, d645069 , to make this PR more compatible with other versions of aiida-core?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, good point, we can do that, I don't think there should be something incompatible in doing so. We just need to open an issue so we don't forget to make the code cleaner. Or actually, we can even think to override within this PR the get_pymathen_structure?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you could in principle create the sites as PeriodicSite

sites = [
            PeriodicSite(
                species=site.species,
                coords=site.coords,
                lattice=Lattice(self.hubbard_structure.cell, pbc=self.hubbard_structure.pbc),
                coords_are_cartesian=True
            ) for site in self.hubbard_structure.get_pymatgen().sites
        ]

And then instead of using the pymat StructureData , just call the function from pymatgen, that it accepts a list of periodic sites ?

get_neighbor_list(r: float, sites: Sequence[PeriodicSite] | None = None, numerical_tol: float = 1e-08, exclude_self: bool = True)→ tuple[np.ndarray, ...][source]

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants