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

Periodicity: StructureData #907

Merged

Conversation

AndresOrtegaGuerrero
Copy link
Contributor

This PR implements changes in the base, and relax workchain so it can create modifications when dealing with 1D or 2D StructureData. This implementation addresses the issues highlighted in #896

@AndresOrtegaGuerrero AndresOrtegaGuerrero mentioned this pull request Apr 13, 2023
@AndresOrtegaGuerrero
Copy link
Contributor Author

@mbercx let me know what you think

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.

Sorry @AndresOrtegaGuerrero! Had to look into the cell_dofree options a bit, since I'm not super familiar, I've left some questions. :)

As for the implementation, maybe it's a matter of preference, but I'm more of a fan of using dictionaries to deal with "cases" instead of if-elif-else statements:

            pbc_cell_dofree_map = {
                (True, True, True): 'all',
                (True, False, False): 'x',
                (False, True, False): 'y',
                (False, False, True): 'z',
                (True, True, False): '2Dxy',
            }
            try:
                base.pw.parameters['CELL']['cell_dofree'] = pbc_cell_dofree_map[structure.pbc]
            except KeyError:
                raise ValueError(
                    f'Structures with periodic boundary conditions `{structure.pbc}` are not supported.'
                )

@@ -180,6 +180,10 @@ def get_builder_from_protocol(
parameters['SYSTEM']['ecutwfc'] = cutoff_wfc
parameters['SYSTEM']['ecutrho'] = cutoff_rho

#If 2D periodicity set assume_isolate
if (True, True, False) == structure.pbc:
parameters['SYSTEM']['assume_isolated'] = '2D'
Copy link
Member

Choose a reason for hiding this comment

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

I see, so this correction can only be used when the first two lattice vectors are entirely in the x-y plane. So running with (False, True, True) should not be allowed for this reason? Should the third lattice vector be perpendicular to this plane?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@mbercx Thank you for the question. So I included this option given the suggestion from @demiranda-gabriel , Usually If i were to do a slab simulation i would omit this part and just be sure that my structure has enough vacuum in the perpendicular direction of the slab. According to the manual this option is only for xy-plane. Though I am pro removing this line.

Copy link
Member

Choose a reason for hiding this comment

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

Hmm, since the get_builder_from_protocol method is meant to provide a set of sensible defaults, it might be reasonable to keep it since the user can always override it. One sentence in the QE documentation still worries me:

Total energy, forces and stresses are computed in a two-dimensional framework.

@demiranda-gabriel So QE would not consider any z-components of the forces? Wouldn't we want to include these when e.g. optimizing the surface layers?

Copy link
Member

Choose a reason for hiding this comment

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

Nevermind, it clearly does. The initial forces and stresses:

     Forces acting on atoms (cartesian axes, Ry/au):

     atom    1 type  1   force =     0.00000000    0.00000000   -0.01049299
     atom    2 type  1   force =     0.00000000    0.00000000    0.01049299
          total   stress  (Ry/bohr**3)                   (kbar)     P=       -4.00
  -0.00003190  -0.00000000   0.00000000           -4.69       -0.00        0.00
  -0.00000000  -0.00003190   0.00000000           -0.00       -4.69        0.00
   0.00000000   0.00000000  -0.00001767            0.00        0.00       -2.60

are reduced to

     Forces acting on atoms (cartesian axes, Ry/au):

     atom    1 type  1   force =     0.00000000    0.00000000    0.00047838
     atom    2 type  1   force =     0.00000000    0.00000000   -0.00047838
          total   stress  (Ry/bohr**3)                   (kbar)     P=        0.07
   0.00000070   0.00000000   0.00000000            0.10        0.00        0.00
   0.00000000   0.00000070   0.00000000            0.00        0.10        0.00
   0.00000000   0.00000000   0.00000006            0.00        0.00        0.01

Input/output cell:

In [1]: relax_wc = load_node(16355)

In [2]: relax_wc.inputs.structure.cell
Out[2]: 
[[3.6086769266, 0.0, 0.0],
 [-1.8043384633, 3.1252058924, 0.0],
 [0.0, 0.0, 21.3114930844]]

In [3]: relax_wc.outputs.output_structure.cell
Out[3]: 
[[3.5615977430325, 0.0, 0.0],
 [-1.7807988715162, 3.0844341234181, 0.0],
 [0.0, 0.0, 21.311492898049]]

if x_dim and not (y_dim or z_dim):
base.pw.parameters['CELL']['cell_dofree'] = 'x'
elif y_dim and not (x_dim or z_dim):
base.pw.parameters['CELL']['cell_dofree'] = 'y'
Copy link
Member

@mbercx mbercx Apr 14, 2023

Choose a reason for hiding this comment

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

Maybe you are more of an expert, but I'm wondering: assume you have a structure with pbc == [False, True, False], and the following CELL_PARAMETERS:

CELL_PARAMETERS angstrom
      13.866974650       0.0000000000       0.0000000000
      1.9334873250       3.3488982827       0.0000000000
      1.9334873250       1.1162994276       13.157371580

I’m wondering if I set cell_dofree to 'y', will only the second component of the second lattice parameter (i.e. 3.3488982827) be allowed to change? If only the second component is allowed to change, also the direction of the lattice vector will change, which could deform the chain structure unless all atoms are on the second lattice vector, no?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

In this scenario , and according to the manual 'y' : only the y component of axis 2 (v2_y) is moved given the cell parameters you are giving there is an angle in that direction so perhaps your pbe = True, False, True is not what you want

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@mbercx Perhaps if there are so many options for cell_dofree the best for QEApp, is just to define the StructureData.pbc (via widget) and just use the logic from this PR externally (not in the RelaxWorkchain) to override the QEworkchain. Since for the QEApp we might just need to consider only a few cases that are 1D (along one direction x,y,z) , 2D (in xy) and 3D. In anycase I know that create_kpoints_from_distance will take into account StructureData.pbc to generate the mesh. Let me know your opinion I can close the PR

Copy link
Member

Choose a reason for hiding this comment

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

In this scenario , and according to the manual 'y' : only the y component of axis 2 (v2_y) is moved given the cell parametersyou are giving there is an angle in that direction so perhaps yourpbe = [True, False, True]` is not what you want

Yeah, it seems you typically would make sure the Xth lattice vector only has an X component. E.g.

CELL_PARAMETERS angstrom
-24.711115283      0.0000            0.00000
0.00000           -24.808651431      0.00000
0.00000            0.000000          4.165100096

Would be a reasonable cell with cell_dofree = 'z'. So a user has to be careful, but we can also help by adding validation on the PwCalculation class (but I won't subject you to that ^^).

@mbercx Perhaps if there are so many options for cell_dofree the best for QEApp, is just to define the StructureData.pbc (via widget) and just use the logic from this PR externally (not in the RelaxWorkchain) to override the QEworkchain. Since for the QEApp we might just need to consider only a few cases that are 1D (along one direction x,y,z) , 2D (in xy) and 3D. In anycase I know that create_kpoints_from_distance will take into account StructureData.pbc to generate the mesh. Let me know your opinion I can close the PR

I understand this will probably be quicker for your use case, and it possible avoids creating bugs in the main aiida-quantumespresso plugin that other users could suffer from. However, long term I think it's definitely something that should be implemented in the aiida-quantumespresso plugin, and since it relies on the pbc it doesn't affect any user running 3D structures at all. So I'd be inclined to push through and get it in here, with a rapid release promised by yours truly.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@mbercx Ok, I will do the changes you suggested, I am pro removing the parameters['SYSTEM']['assume_isolated'] = '2D'

Copy link
Member

Choose a reason for hiding this comment

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

@AndresOrtegaGuerrero if it's ok for you, I was already doing the changes and adding a test, but got distracted by aiidateam/aiida-core#5967 when I actually tried to import a 2D structure. 😅

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@mbercx no problem! , thank you

Copy link
Member

Choose a reason for hiding this comment

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

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 a lot @AndresOrtegaGuerrero! This LGMT 🚀

@mbercx mbercx merged commit 5687800 into aiidateam:main Apr 17, 2023
mbercx added a commit that referenced this pull request Apr 17, 2023
Currently none of the `get_builder_from_protocol()` methods consider the
periodic boundary conditions (`pbc`) of the input `StructureData`. Here
we add support for several cases to the `PwBaseWorkChain` and
`PwRelaxWorkChain`:

* 2D structures in the x-y plane with periodic boundary conditions for
  the first two lattice vectors (`pbc == [True, True, False]`). For this
  case, the `PwBaseWorkChain` protocol will set `SYSTEM.assume_isolated`
  to `2D`, and the `PwRelaxWorkChain` protocol will use
  `CELL.do_free = '2Dxy'` in case the `RelaxType` is `CELL` or
  `POSITIONS_CELL`.
* 1D structures defined on an arbitrary lattice vector will use the
  corresponding `CELL.do_free` setting (`'x'`, `'y'`, `'z'`) if
  `RelaxType` is `CELL` or `POSITIONS_CELL`.

So far, no validation is done to see if the structure has the correct
`cell` specified. This could potentially be added to the `PwCalculation`
validation.

Co-authored-by: Marnik Bercx <mbercx@gmail.com>
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.

2 participants