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

adapting BandsWorkChain for HubbardStructureData #998

Merged

Conversation

AndresOrtegaGuerrero
Copy link
Collaborator

@AndresOrtegaGuerrero AndresOrtegaGuerrero commented Dec 16, 2023

When using HubbardStructureData with PwBandsStructure the workchain doesnt use the hubbard information.
Since it calls seekpath , the information is lost.

Additionally the HubbardStructureData has a bug for StructureData with periodicity 1d and 2d.

This PR addresses these two issues , #999

@bastonero
Copy link
Collaborator

Thanks @AndresOrtegaGuerrero. I am wondering if then all the methods are working as expected. Could you please provide some tests and/or examples, with some structure? I think there should be already some fixtures for 2D. It is important that the functions work, and in the intended behaviour. Here, it is not clear if your fix is going to work out.

@AndresOrtegaGuerrero
Copy link
Collaborator Author

AndresOrtegaGuerrero commented Dec 23, 2023

Hi @bastonero , The 2D and 1D is already working at the PwBaseWorkChain (like the sampling of k-points) we did these with @mbercx. Regarding the PwBandsWorkChain, this branch manage to fix the bug. Another test is the PdosWorkChain that works without no problems with HubbardStructureData. Is also tested in the AiiDAlab qe app (aiidalab/aiidalab-qe#577)

@AndresOrtegaGuerrero
Copy link
Collaborator Author

I think you meant to add test for the hubbard structure, and for bands workchain 😅

Copy link
Member

@unkcpz unkcpz left a comment

Choose a reason for hiding this comment

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

@AndresOrtegaGuerrero thanks! I suggest to add a get_primitive_structure to HubbardStructureData and call it to extract the bind information from a single HubbardStructureData object. See the comment below.

I am not super familiar what habbard_info contains, maybe it doesn't relate to whether the structure is primitive or not. So the change may not needed.

Comment on lines 114 to 116
ase = self.get_ase()
ase.set_pbc((True, True, True))
temp_structure = StructureData(ase=ase)
Copy link
Member

Choose a reason for hiding this comment

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

I think it should be able to create a temp structure directly with StructureData (https://github.com/aiidateam/aiida-core/blob/06ea130df8854f621e25853af6ac723c37397ed0/src/aiida/orm/nodes/data/structure.py#L1592) instead of through ase structure.

What will happened if the pbc is not 3D? It will raise an error or it will give the wrong results afterwards?

Copy link
Collaborator

Choose a reason for hiding this comment

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

why not tmp_hubbard_structure = hubbard_structure.clone() for instance, and then temp_hubbard_structure.pbc = (...)? At least there is a good chance we are not missing anything.

Copy link
Collaborator

Choose a reason for hiding this comment

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

moreover, this is lacking of a test asserting that the pymatgen module would work fine and produce the expected translation vectors, crucial for intersite parameters. and, very important, we should not lose the kind names and anything while doing structuredata --> ase --> structuredata, so that's why the approach i suggested in the previous comment. when used with only onsite hubbard U it can be fine, but the all point here is indeed to have the V parameters set correctly with translation vectors, so please add a careful test, otherwise this would break the implementation for extended Hubbard parameters

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

if is different that 3D , it will try to use pymatgen, and then it will fail

Copy link
Collaborator Author

@AndresOrtegaGuerrero AndresOrtegaGuerrero Jan 28, 2024

Choose a reason for hiding this comment

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

@unkcpz i cant use set_pbc since , is an StructureData that is saved

Copy link
Member

Choose a reason for hiding this comment

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

Did you try to pass the pbc attribute when create StructureData? I assume you can do:

StructureData(pbc=[True, True, True])

Copy link
Collaborator Author

@AndresOrtegaGuerrero AndresOrtegaGuerrero Jan 31, 2024

Choose a reason for hiding this comment

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

you mean like this temp_structure = StructureData(ase=ase, pbc=[True, True, True]) ? because since the structure data is stored i cant do any modifications

Copy link
Member

Choose a reason for hiding this comment

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

No I mean using as less structure type as possible.
With your implementation, the structure type used are HubbardStructureData -> asc structure -> StructureData -> pymatgen structure. However, the purpose here is to get the sites which are information that should include in the StructureData.

I think we should try our best to not use ase/pymatgen, at least should try to not add more dependent on those. I want to see @mbercx's opinion on it.

Copy link
Member

@mbercx mbercx Jan 31, 2024

Choose a reason for hiding this comment

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

Thanks for the ping @unkcpz.

I see your point that going through ASE to get a StructureData again to then get a pymatgen structure seems very roundabout, but I understand why @AndresOrtegaGuerrero did it.

The issue is that get_pymatgen() would not work since it would return a pymatgen Molecule, and get_pymatgen_structure() fails since that only allows periodic boundary conditions... Ideally this would actually be possible. I don't immediately understand why StructureData doesn't allow it, since a pymatgen Structure can also e.g. have (True, False, False) periodic boundary conditions. I'll open an issue on aiida-core to discuss this, but that'll take time.

Unfortunately, it's also not so easy to get an (unstored) copy of a StructureData directly. We could avoid using ASE in another way though. In the end the only thing the method needs is the sites, so the following should work for structures of any dimension:

from pymatgen.core import PeriodicSite, Lattice

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

And is probably cleaner than the if/else clause with the StructureData -> ASE -> StructureData -> Structure -> sites route for <3D structures.

As a final note, I know @unkcpz wants to avoid using ASE/pymatgen since these can be somewhat problematic dependencies (ASE because there are barely any releases, pymatgen since they often make breaking changes). But in the end I don't see a way to avoid using these tools, nor do I think it is desirable. That would involve reimplementing (and maintaining!) a lot of functionality ourselves.

if isinstance(self.ctx.current_structure, HubbardStructureData):
self.ctx.hubbard_info = self.ctx.current_structure.hubbard
hubbard_structure = HubbardStructureData.from_structure(result['primitive_structure'])
hubbard_structure.hubbard = self.ctx.hubbard_info
Copy link
Member

Choose a reason for hiding this comment

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

I am a bit worried here, the habbard_info is bind with the original structure. But the current_structure will then be a primitive structure.

I think it will be more proper that there will be a get_primitive_structure method for HubbardStructureData and hubbard_info is always attached with the corresponded structure.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Just passing by: be very aware that getting a primitive structure would mean to also transfer the hubbard parameters CORRECTLY, e.g. in case of extended Hubbard parameters. This is still not implemented in e.g. HubbardUtils, for which there is a method that is only tested for supercells. And indeed, this way would most certainly set WRONG hubbard parameters for the primitive structure, where the order of atoms might change, as well as the number, and considering different lattice parameters this is a recepy for disaster. As such, is it REALLY necessary to have a primite structure? If not, then it would be wiser to just skip this if uncertain on what it is doing.

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 wouldnt change the structure, since this function is just to pass the hubbard parameters, instead i proposed that 1D and 2D , are treated as 3D , just to pass the parameters to the HubbardStructureData,

Copy link
Collaborator

Choose a reason for hiding this comment

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

Unfortunately, the primitive_structure can be different from the current_structure. The issue is that the hubbard parameters are defined on top of the current_structure, and you cannot pass such parameters as they are to the primitive_structure.

Simple example where this implementation wouldn't work. Say you have a current_structure with 4 atoms, namely Fe, Fe, O, O. Let them be ordered in the current_structure.sites, so that current_structure.sites[0] corresponds to the first Fe and so on. Now, the hubbard datastructure would save the hubbard parameters using the indices to define the parameters, hence U(Fe)=5.0 will be defined by the tuple (0,'3d',0,'3d',(0,0,0), 5.0), for instance. Let's assume that you also have a U parameter on all the other atoms, meaning also the oxygen would have e.g. (2,'2p',2,'2p',(0,0,0), 8.0). If now your primitive_structure, being primitive, is smaller than the current_structure, meaning e.g. it has only Fe,O in its sites attribute (because the other two atoms of the current_structure are symmetry equivalent), than you see that you cannot pass anymore the hubbard information, as you don't have anymore an index 2 in primitive_structure.sites. Moreover, even if you still have the same number of sites, it's not granted that the order of the atoms (in the sites) will be the same as in the current_structure, and the hubbard parameters would be wrongly assigned in this naive way. Not to speak when one has an intersite interaction, e.g. (0,'3d',2,'2p',(1,1,1), 8.0), where (1,1,1) is the translation vector defined with respect to the current_structure cell and not to the primitive_structure one, hence completely breaking the intersite information.

Hope this clarifies that it's not as trivial as it seems to transfer such parameters. I suggest you to look into the implementation, and e.g. have a look to aiida_quantumespresso.utils.hubbard and in particular to the HubbardUtils.get_hubbard_for_supercell method, to see how the reindexing can be far from trivial.

This said, probably it would be just easier to just skip the use of a primitive_structure, which should not be strictly necessary for having an electronic band structure.

Copy link
Member

Choose a reason for hiding this comment

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

From the docstring of seekpath_structure_analysis, it says the calcfuntion will "Primitivize the structure with SeeKpath". Can you test on a conventional cell with the code you implement and see if the hubbard structure is changed and still valid to conform with the hubbard_info.
As I know, the band structure should be conduct on the primitive cell to avoid the band folded in BZ, and that is the reason the result['primitive_structure'] is used.

Copy link
Collaborator Author

@AndresOrtegaGuerrero AndresOrtegaGuerrero Jan 31, 2024

Choose a reason for hiding this comment

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

@unkcpz Indeed when using seekpath and generating a primitive_structure , we CAN'T use hubbard_info. Is also quite tricky to generalize it since the HubbardStructureData can have hubbard_parameters for DFT+U, DFT+U+V (onsites and intersites parameters). It requires a method to generalize the manifold to apply it on a primitive_structure

Copy link
Member

Choose a reason for hiding this comment

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

Then can you make it except when the scenario happened? If I understand correctly, with your implementation, if I run the band on the supercell with hubbard turned on, it will raise in some place that the calculation will fail. Can you try to add a validation and notify user about it.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@unkcpz I added a logic that generalize the hubbard_info , to be transformed to the primitive_structure. @bastonero could you give me your opinion , I will try to make a test with a conv unit cell , asap

inputs = {
'reference_distance': self.inputs.get('bands_kpoints_distance', None),
'metadata': {
'call_link_label': 'seekpath'
}
}
result = seekpath_structure_analysis(self.ctx.current_structure, **inputs)
self.ctx.current_structure = result['primitive_structure']

# If the input structure was a HubbardStructureData
Copy link
Member

Choose a reason for hiding this comment

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

I haven't check the logic here, but instead of adding this to the PwBandsWorkChain, wouldn't it make more sense to move this into the seekpath_structure_analysis?

Copy link
Collaborator

Choose a reason for hiding this comment

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

I would opt for placing this in the HubbardUtils, so that we can reuse the logic elsewhere.
The logic for U seems fine, for V unfortunately it's more complex. Moreover, the initialize_intersites_hubbard is not really meant for this, but one should use the append_hubbard_parameters one. The difficulties are the following:

  1. The number of atoms in the unitcell are greater than the primitive, hence the indices won't be the same
  2. The V interaction is specific to a pair of atoms, that's why the initialize_intersites_hubbard is not suited for this
  3. One needs to find the translation vectors and the corresponding atoms by symmetry in the primitive cell

A good proxy for testing is to use the HubbardUtils.get_hubbard_for_supercell once you have the new paramters in the primitive. The routine should take the unitcell as input, and the returned hubbard_structure should have the same hubbard parameters as the original.

Please, have a look into HubbardUtils.get_hubbard_for_supercell to understand the logic. Probably it would be sufficient to copy/paste this routine and adapt the logic for the primitive cell. The current method wouldn't work as the primitive cell is smaller, hence the translation vectors would be fractional, and the logic doesn't account for this. One would probably need to use the inverse cell of the primitive, so to "restore" integer translations.

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 tried the logic, it worked , should i move it to seekpath_structure_analysis calcfunction? (I am ok with :D)

Copy link
Collaborator

Choose a reason for hiding this comment

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

Thanks a lot Andres!
What do you mean it worked? Does this logic also work for U+V, is that what you mean?

@@ -61,6 +61,19 @@ def test_from_structure(generate_structure, generate_hubbard):
assert len(hubbard_structure.kinds) == 2


@pytest.mark.usefixtures('aiida_profile')
@pytest.mark.parametrize('structure_name', ['2D-xy-arsenic', '1D-x-carbon', '1D-y-carbon', '1D-z-carbon'])
Copy link
Member

Choose a reason for hiding this comment

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

This test doesn't seem to test the append_hubbard_parameter method, which is the one you are actually adapting?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

It should, since the test is to check that it doesnt fail with pbc different that 3D

@mbercx
Copy link
Member

mbercx commented Jan 31, 2024

Thanks @AndresOrtegaGuerrero! One more general note here is that you are making two quite different changes in one PR:

  1. Provide support for low-dimensional structures to HubbardStructureData.
  2. Provide support for HubbardStructureData to the PwBandsWorkChain.

No need to open another PR, but I'd split these changes up into two commits that describe each. Happy to do this once the PR is happy to merge, or discuss how to do this over Zoom.

@AndresOrtegaGuerrero
Copy link
Collaborator Author

Thanks @AndresOrtegaGuerrero! One more general note here is that you are making two quite different changes in one PR:

  1. Provide support for low-dimensional structures to HubbardStructureData.
  2. Provide support for HubbardStructureData to the PwBandsWorkChain.

No need to open another PR, but I'd split these changes up into two commits that describe each. Happy to do this once the PR is happy to merge, or discuss how to do this over Zoom.

Thank you @mbercx sorry about the two changes, in my head they were quite small that since it was an aiidalab_qe issue i thought it was ok in one. It wont happen again :D

@mbercx mbercx force-pushed the fix/hubbard_bands_periodicity branch 2 times, most recently from c5a3556 to b51518a Compare February 2, 2024 10:14
@AndresOrtegaGuerrero AndresOrtegaGuerrero force-pushed the fix/hubbard_bands_periodicity branch from b51518a to b29b561 Compare February 2, 2024 11:51
@bastonero
Copy link
Collaborator

The pbc should be fixed by this PR aiidateam/aiida-core#6281 (comment)

@AndresOrtegaGuerrero
Copy link
Collaborator Author

@bastonero , @mbercx proposed a solution here without the need of using generate_structure

Copy link
Member

@mbercx mbercx left a comment

Choose a reason for hiding this comment

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

Thanks @AndresOrtegaGuerrero. As @bastonero already mentioned, I think for intersite Hubbard parameters it's not so trivial to correctly preserve the parameters for the primitive structure provided by SeeK-path (see comments).

To still allow the user to find the standardized path in the intersite V case, we should see if we can't fully preserve the structure and have SeeK-path return the corresponding path. You'd have to see if this is possible, maybe touch base with @giovannipizzi on the topic.

return _generate_hubbardstructure_conv


# pylint: disable=W0621
Copy link
Member

Choose a reason for hiding this comment

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

I think it's clearer to use the string identifier for pylint disables.

Suggested change
# pylint: disable=W0621
# pylint: disable=redefined-outer-name


# pylint: disable=W0621
@pytest.mark.usefixtures('aiida_profile')
def test_seekpath_analysis(generate_hubbardstructure_conv):
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
def test_seekpath_analysis(generate_hubbardstructure_conv):
def test_hubbard(generate_hubbardstructure_conv):

return result


def update_structure_with_hubbard(structure, result):
Copy link
Member

Choose a reason for hiding this comment

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

I've had a closer look, and I'm wondering if the current approach is correct. As @bastonero mentions here, primitivizing a HubbardStructureData and correctly passing the parameters can be quite challenging.

After a bit of testing, it seems SeeKpath does respect the kinds when primitivizing (makes sense, also e.g. for magnetic runs this would be critical):

from aiida import orm, load_profile
from ase import build
from aiida_quantumespresso.calculations.functions.create_magnetic_configuration import create_magnetic_configuration
from aiida.tools.data.array.kpoints import get_explicit_kpoints_path

load_profile()

conv = orm.StructureData(ase=build.bulk('Co', crystalstructure='fcc', a=5.43, cubic=True))
new_struc = create_magnetic_configuration(conv, [0, 0, 1, 1], metadata={'store_provenance': False})['structure']
get_explicit_kpoints_path(new_struc)['primitive_structure'].kinds

So for onsite parameters (U) your approach should work fine, as defining different U parameters for two sites would require different kinds as well, and I think your approach works as intended.

For intersite (V) parameters, I'm quite sure the current approach is incomplete. Admittedly, I'm not very experienced with these calculations, so having @bastonero's input would be great. As a basic test, I was curious what would happen if I took an example structure which is already primitive and use the new seekpath_structure_analysis. Here is the Hubbard card before:

HUBBARD	ortho-atomic
 V	Ni-3d	Ni-3d	1	1	9.6168
 V	Ni-3d	N-2p	1	3	1.034
 V	Ni-3d	N-2p	1	69	1.034
 V	Ni-3d	N-2p	1	51	1.034
 V	Ni-3d	N-2p	1	78	1.034
 V	Ni-3d	Ni-3d	2	2	9.5156
 V	Ni-3d	N-2p	2	3	1.1805
 V	Ni-3d	N-2p	2	45	1.1805

and after:

HUBBARD	ortho-atomic
 V	Ni-3d	Ni-3d	1	1	9.6168
 V	Ni-3d	N-2p	1	3	1.034
 V	Ni-3d	Ni-3d	1	1	9.5156
 V	Ni-3d	N-2p	1	3	1.1805

Which now seems to define a different V interaction for the exact same pair of sites, so something seems wrong here. I think the issue is that the kinds are no longer enough to capture the pair of sites to which apply the V. You'd at least have to use the atom indices, and those aren't even enough because you can have a different intersite interaction with the same site in different neighboring cells (defined via the translations). So this is a whole other ball game.

So maybe for now it'd be best to only support onsite interactions?

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 guess for V, we will need to check distancee between index of atoms , and them compare and find the distance of the Hubbard(inital) to find which atoms in the primitive have the same distance to assign the new V, and in this case we shouldnt use initialize_intersites_hubbard, since the Hubbard you provide in your example considers a different V depending on the distance

Copy link
Member

Choose a reason for hiding this comment

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

I guess for V, we will need to check distancee between index of atoms , and them compare and find the distance of the Hubbard(inital) to find which atoms in the primitive have the same distance to assign the new V

Would the distance between sites be sufficient? Can't there be two pairs of atoms that have the same distance but a different intersite V? In that case we'd have to check that the primitivization respects each V pair, not just within the unit cell using the atom index, but also the neighbors via the translations.

and in this case we shouldnt use initialize_intersites_hubbard, since the Hubbard you provide in your example considers a different V depending on the distance

Indeed, append_hubbard_parameter would allow you to specify the atom/site indices and translation.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Indeed append_hubbard_parameters is suited for that. You can avoid to specify the translation, and the routine will find the image of the V atom closest to the U atom, i.e. its neighbour. I bet this is sufficient, so one doesn't have to specify the translation vector, which would be messy.
Also note: the indices should be given within primitive cell, so one should get from the seekpath analysis the mapping between atom indices from unit cell to primitive cell, so that for each HubbardParameter you can map the indices correspondingly.
Nevertheless, a better logic that also is able to specify the exact translation vector would be better. As I already suggested, one could follow the logic of the routine I implemented in HubbardUtils.


def _generate_hubbardstructure_conv():
cell = [[7.96416, 0.0, 4.87664153e-16], [-4.87664153e-16, 7.96416, 4.87664153e-16], [0.0, 0.0, 7.96416]]
sites = [['Li', 'Li', (2.98656, 0.99552, 6.96864)], ['Li', 'Li', (6.96864, 6.96864, 4.9776)],
Copy link
Member

Choose a reason for hiding this comment

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

As already mentioned, this structure seems excessively large for this test. Besides this, I'd like to see more tests that the primitive structure has the correct Hubbard parameters. If there is a bug here, other users might start mindlessly using the seekpath_structure_analysis calculation function and not notice that the resulting HubbardStructureData has incorrect parameters.

Currently, you're only really testing the case where all the sites of the same element have the same kind name. I'd still want to make sure the current approach works for the case where e.g. there are two Co kinds with different U values. You wouldn't need such a large cell for this, a simple conventional FCC Co structure would to the trick.

assert isinstance(
result['primitive_structure'], HubbardStructureData
), 'Primitive structure should be a HubbardStructureData'
assert conventional_parameters != primitive_parameters, 'Primitive parameters should be different'
Copy link
Member

Choose a reason for hiding this comment

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

They should only really be different in case the structure is different. I'd also like to see that the parameters are the same when the structure is the same. In many cases I'm running an already primitivized structure, so running SeeKpath would not change the structure (this isn't always true, there can be cell rotations).

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 guess we can check the logic, if the structure is different from the primitive seekpath is returning ?

assert conventional_parameters != primitive_parameters, 'Primitive parameters should be different'
assert len(primitive_parameters) == len(
conventional_parameters
), 'Primitive parameters should have the same length as conventional parameters'
Copy link
Member

Choose a reason for hiding this comment

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

This test would fail for my intersite Ni2N example in the comment above.

assert all(
conv_param.atom_manifold == prim_param.atom_manifold
for conv_param, prim_param in zip(conventional_parameters, primitive_parameters)
), 'Atom manifold should match in primitive'
Copy link
Member

Choose a reason for hiding this comment

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

Again, I think this is only really testing the onsite case. To check if the seekpath_structure_analysis works in all cases, we'd have to be much more rigorous.

@mbercx mbercx force-pushed the fix/hubbard_bands_periodicity branch 3 times, most recently from 3389e75 to 40de621 Compare February 12, 2024 17:41
@mbercx
Copy link
Member

mbercx commented Feb 12, 2024

@AndresOrtegaGuerrero @bastonero I'll move the seekpath-related changes to a different PR so we can focus on the <3D structures here.

I've added some validation on the site indices for append_hubbard_parameter in f8d6933.

I've also extended the tests of the method to check more cases and verify that adapting it for <3D structures doesn't break anything. ^^

@mbercx mbercx force-pushed the fix/hubbard_bands_periodicity branch from 40de621 to 72741f6 Compare February 12, 2024 17:50
@mbercx
Copy link
Member

mbercx commented Feb 12, 2024

@bastonero do you still want to have a look?

@bastonero
Copy link
Collaborator

@bastonero do you still want to have a look?

Let's have a quick one

Comment on lines 114 to 116
ase = self.get_ase()
ase.set_pbc((True, True, True))
temp_structure = StructureData(ase=ase)
Copy link
Collaborator

Choose a reason for hiding this comment

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

why not tmp_hubbard_structure = hubbard_structure.clone() for instance, and then temp_hubbard_structure.pbc = (...)? At least there is a good chance we are not missing anything.

tests/data/test_hubbard_structure.py Show resolved Hide resolved
Comment on lines 114 to 116
ase = self.get_ase()
ase.set_pbc((True, True, True))
temp_structure = StructureData(ase=ase)
Copy link
Collaborator

Choose a reason for hiding this comment

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

moreover, this is lacking of a test asserting that the pymatgen module would work fine and produce the expected translation vectors, crucial for intersite parameters. and, very important, we should not lose the kind names and anything while doing structuredata --> ase --> structuredata, so that's why the approach i suggested in the previous comment. when used with only onsite hubbard U it can be fine, but the all point here is indeed to have the V parameters set correctly with translation vectors, so please add a careful test, otherwise this would break the implementation for extended Hubbard parameters

if isinstance(self.ctx.current_structure, HubbardStructureData):
self.ctx.hubbard_info = self.ctx.current_structure.hubbard
hubbard_structure = HubbardStructureData.from_structure(result['primitive_structure'])
hubbard_structure.hubbard = self.ctx.hubbard_info
Copy link
Collaborator

Choose a reason for hiding this comment

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

Unfortunately, the primitive_structure can be different from the current_structure. The issue is that the hubbard parameters are defined on top of the current_structure, and you cannot pass such parameters as they are to the primitive_structure.

Simple example where this implementation wouldn't work. Say you have a current_structure with 4 atoms, namely Fe, Fe, O, O. Let them be ordered in the current_structure.sites, so that current_structure.sites[0] corresponds to the first Fe and so on. Now, the hubbard datastructure would save the hubbard parameters using the indices to define the parameters, hence U(Fe)=5.0 will be defined by the tuple (0,'3d',0,'3d',(0,0,0), 5.0), for instance. Let's assume that you also have a U parameter on all the other atoms, meaning also the oxygen would have e.g. (2,'2p',2,'2p',(0,0,0), 8.0). If now your primitive_structure, being primitive, is smaller than the current_structure, meaning e.g. it has only Fe,O in its sites attribute (because the other two atoms of the current_structure are symmetry equivalent), than you see that you cannot pass anymore the hubbard information, as you don't have anymore an index 2 in primitive_structure.sites. Moreover, even if you still have the same number of sites, it's not granted that the order of the atoms (in the sites) will be the same as in the current_structure, and the hubbard parameters would be wrongly assigned in this naive way. Not to speak when one has an intersite interaction, e.g. (0,'3d',2,'2p',(1,1,1), 8.0), where (1,1,1) is the translation vector defined with respect to the current_structure cell and not to the primitive_structure one, hence completely breaking the intersite information.

Hope this clarifies that it's not as trivial as it seems to transfer such parameters. I suggest you to look into the implementation, and e.g. have a look to aiida_quantumespresso.utils.hubbard and in particular to the HubbardUtils.get_hubbard_for_supercell method, to see how the reindexing can be far from trivial.

This said, probably it would be just easier to just skip the use of a primitive_structure, which should not be strictly necessary for having an electronic band structure.

inputs = {
'reference_distance': self.inputs.get('bands_kpoints_distance', None),
'metadata': {
'call_link_label': 'seekpath'
}
}
result = seekpath_structure_analysis(self.ctx.current_structure, **inputs)
self.ctx.current_structure = result['primitive_structure']

# If the input structure was a HubbardStructureData
Copy link
Collaborator

Choose a reason for hiding this comment

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

I would opt for placing this in the HubbardUtils, so that we can reuse the logic elsewhere.
The logic for U seems fine, for V unfortunately it's more complex. Moreover, the initialize_intersites_hubbard is not really meant for this, but one should use the append_hubbard_parameters one. The difficulties are the following:

  1. The number of atoms in the unitcell are greater than the primitive, hence the indices won't be the same
  2. The V interaction is specific to a pair of atoms, that's why the initialize_intersites_hubbard is not suited for this
  3. One needs to find the translation vectors and the corresponding atoms by symmetry in the primitive cell

A good proxy for testing is to use the HubbardUtils.get_hubbard_for_supercell once you have the new paramters in the primitive. The routine should take the unitcell as input, and the returned hubbard_structure should have the same hubbard parameters as the original.

Please, have a look into HubbardUtils.get_hubbard_for_supercell to understand the logic. Probably it would be sufficient to copy/paste this routine and adapt the logic for the primitive cell. The current method wouldn't work as the primitive cell is smaller, hence the translation vectors would be fractional, and the logic doesn't account for this. One would probably need to use the inverse cell of the primitive, so to "restore" integer translations.

tests/calculations/functions/test_seekpath_analysis.py Outdated Show resolved Hide resolved
tests/data/test_hubbard_structure.py Outdated Show resolved Hide resolved
return result


def update_structure_with_hubbard(structure, result):
Copy link
Collaborator

Choose a reason for hiding this comment

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

Indeed append_hubbard_parameters is suited for that. You can avoid to specify the translation, and the routine will find the image of the V atom closest to the U atom, i.e. its neighbour. I bet this is sufficient, so one doesn't have to specify the translation vector, which would be messy.
Also note: the indices should be given within primitive cell, so one should get from the seekpath analysis the mapping between atom indices from unit cell to primitive cell, so that for each HubbardParameter you can map the indices correspondingly.
Nevertheless, a better logic that also is able to specify the exact translation vector would be better. As I already suggested, one could follow the logic of the routine I implemented in HubbardUtils.

src/aiida_quantumespresso/data/hubbard_structure.py Outdated Show resolved Hide resolved
@bastonero
Copy link
Collaborator

Lol, I had a huge review with all the comments that nobody was seeing. That's why nobody was replying to those comments

@mbercx
Copy link
Member

mbercx commented Feb 12, 2024

Lol, I had a huge review with all the comments that nobody was seeing. That's why nobody was replying to those comments

Haha, I read "let's have a quick one", then saw the massive review and almost fell out of my chair. 🤣

@mbercx
Copy link
Member

mbercx commented Feb 12, 2024

Now working on #1007, will consider your comments there, but we already agreed that we'd stick to the onsite case for now.

Copy link
Collaborator

@bastonero bastonero left a comment

Choose a reason for hiding this comment

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

Thanks @mbercx , I do approve 🚀

@bastonero
Copy link
Collaborator

@mbercx it seems you have to self approve yourself :D

@mbercx
Copy link
Member

mbercx commented Feb 12, 2024

@mbercx it seems you have to self approve yourself :D

Oh right, I requested changes after all 😅

@mbercx mbercx self-requested a review February 12, 2024 21:48
@mbercx
Copy link
Member

mbercx commented Feb 12, 2024

Haha, seems we also need @unkcpz's approval now 😂

mbercx and others added 2 commits February 13, 2024 11:57
The `append_hubbard_parameter()` method of the `HubbardStructureData` class currently
does not verify that the `atom_index` and `neighbour_index` values correspond to valid
site indices in the structure. This means that the method will:

- fail in case `translation` is set to `None` (as it attempts to access these site indices to obtain
  the distance and image of the nearest neighbour using periodic boundary conditions)
- succeed in case the translation is specified, but fail once the code tries to write Hubbard card
  during the `prepare_for_submission` method of the `PwCalculation`.

To have more consistent behaviour and avoid exceptions during the `upload` step of the
submission process, we check here that the `atom_index` and `neighbour_index` are valid and raise
a `ValueError` if they are not.
The `HubbardStructureData.append_hubbard_parameter()` method currently fails for
structures that don't have 3-dimensional periodic boundary conditions. This is caused
by a restriction in the StructureData.get_pymatgen_structure()` method, which raises a
`ValueError` if the periodic boundary conditions are not (True, True, True). This has
been fixed in a recent commit on `aiida-core`:

aiidateam/aiida-core@adcce4b

But this will take some time to be released and would restrict our compatibility to
`aiida-core>=2.6.0`.

Here we directly construct a list of `PeriodicSite` objects from the `sites` list of the
object obtained from `StructureData.get_pymatgen()`, which can be both either a
pymatgen `Structure` or `Molecule`, depending on the dimensionality of the structure.

Co-authored-by: AndresOrtegaGuerrero <[email protected]>
@mbercx mbercx force-pushed the fix/hubbard_bands_periodicity branch from 4b5d9f2 to e190376 Compare February 13, 2024 10:58
Copy link
Member

@unkcpz unkcpz left a comment

Choose a reason for hiding this comment

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

Looks cool, already out of my scope! Thanks everyone, approved.

@mbercx mbercx merged commit d645069 into aiidateam:main Feb 13, 2024
7 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants