This is an Issue to track the development of the CLI functionality for this package. When development on this feature begins, will create a new branch for this and link this Issue.

Would aim to be able to do the following tasks via CLI:

• Structure generation
• Parsing to create csv of outputs with energies vs distortion, RMS etc (Move the parsing_script.sh functionality into this)
• Generate energies vs distortion plots

Hi I got this error with snb-regenerate.

kanta@ln04:~/work/SrTiO3_Al/defect/vac/vac_Sr> snb-regenerate -v

vac_1_Sr
Parsing vac_1_Sr_-2...
vac_1_Sr_-2: Energy difference between minimum, found with Rattled bond distortion, and unperturbed: -0.73 eV.
Energy lowering distortion found for vac_1_Sr with charge -2. Adding to low_energy_defects dictionary.
Parsing vac_1_Sr_-1...
No energy lowering distortion with energy difference greater than min_e_diff = 0.05 eV found for vac_1_Sr with charge -1.
Parsing vac_1_Sr_0...
Bond_Distortion_-60.0% not fully relaxed
vac_1_Sr_0: Energy difference between minimum, found with -0.5 bond distortion, and unperturbed: -0.14 eV.
Comparing structures to specified ref_structure (Sr26 Ti27 O81)...
New (according to structure matching) low-energy distorted structure found for vac_1_Sr_0, adding to low_energy_defects['vac_1_Sr'] list.

Comparing and pruning defect structures across charge states...
io.py:419: UserWarning: /mnt/lustre/a2fs-work1/work/e685/e685/kanta/SrTiO3_Al/defect/vac/vac_Sr/vac_1_Sr_0/Bond_Distortion_-60.0%/CONTCAR file doesn't exist, storing as 'Not converged'. Check path & relaxation
Traceback (most recent call last):
  File "/work/e685/e685/kanta/.local/bin/snb-regenerate", line 8, in <module>
    sys.exit(regenerate())
  File "/work/e685/e685/kanta/.local/lib/python3.8/site-packages/click/core.py", line 1130, in __call__
    return self.main(*args, **kwargs)
  File "/work/e685/e685/kanta/.local/lib/python3.8/site-packages/click/core.py", line 1055, in main
    rv = self.invoke(ctx)
  File "/work/e685/e685/kanta/.local/lib/python3.8/site-packages/click/core.py", line 1404, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/work/e685/e685/kanta/.local/lib/python3.8/site-packages/click/core.py", line 760, in invoke
    return __callback(*args, **kwargs)
  File "/work/e685/e685/kanta/.local/lib/python3.8/site-packages/shakenbreak/cli.py", line 1393, in regenerate
    _ = energy_lowering_distortions.get_energy_lowering_distortions(
  File "/work/e685/e685/kanta/.local/lib/python3.8/site-packages/shakenbreak/energy_lowering_distortions.py", line 621, in get_energy_lowering_distortions
    low_energy_defects = _prune_dict_across_charges(
  File "/work/e685/e685/kanta/.local/lib/python3.8/site-packages/shakenbreak/energy_lowering_distortions.py", line 277, in _prune_dict_across_charges
    comparison_results = compare_struct_to_distortions(
  File "/work/e685/e685/kanta/.local/lib/python3.8/site-packages/shakenbreak/energy_lowering_distortions.py", line 734, in compare_struct_to_distortions
    sorted_distorted_df = matching_sub_df[
TypeError: sort_values() got an unexpected keyword argument 'key'

I am following the doped example workbook as well as the SNB examples page to generate a defect via doped, export to json then reimport using loadfn and inputing into Distortions. The traceback gives:

6 # Create an instance of Distortion class with the defect dictionary and the distortion parameters
----> 8 Dist = Distortions(
      9     defects=doped_defect_dict,  # the defects dictionary
     10     #oxidation_states=oxidation_states,  # explicitly specify oxidation states
     11 )

TypeError: Distortions.__init__() missing 1 required positional argument: 'defects_dict'

I am not sure if this is an issue with doped dumpfn where it's supposed to have a certain section title or if it's because SNB has since been updated as I assume this used to work when the workbook was created.

Any help would be greatly appreciated.

This has been discussed internally, but posting here to keep an online record:

Seán:

Ok data dump incoming. Something I started looking at the other week, was analysing the convergence behaviour of energy and forces vs ionic step in the ShakeNBreak relaxations, to see if there was any correlations that could be leveraged to increase the efficiency of the calcs (e.g. maybe looking at the relative energies after only ~20 ionic steps would be enough to show where energy-lowering distortions are happening, or maybe we could use energy differences rather than forces for EDIFFG to speed things up etc.). Took a while to scrape through the data, but checked through for a few different test cases (V_Cd in CdTe, V_Na in NaBiS2 and Se_i in Se; see attachments).

My conclusions are that

• No real general correlations between energies/forces in the early ionic steps and the final relaxation energy
• No we should def stick with forces rather than energies for EDIFFG
• We can increase our default EDIFFG from -0.01 to -0.025 , which should give a ~30-40% reduction in the number of required ionic steps while still being well-converged for extracting energy differences (i.e. identifying energy lowering distortions)

Note that the decision to increase the default EDIFFG to -0.025 was reversed in #8 after this was found to slightly worsen performance in the case of the tricky V_Se in GeSe, however for high-throughput applications the recommendation will be to use a looser EDIFFG = -0.025 rather than the default EDIFFG = -0.01, as the increase in speed outweighs the small drop in accuracy.

Attached are the files from this analysis, with key screenshots below:

$V_{Cd}$ in CdTe:

Change in energy ($\Delta$ Energy) vs max atomic force per ionic step:

What about changing EDIFFG to energy rather than force-based convergence?

$V_{Na}$ in $NaBiS_2$

$Se_i$ in Se:

E_and_F_vs_Step_SnB_Tests.zip

As mentioned in #2, the new Monte Carlo rattling function allows to include a penalty for interatomic distances below a certain threshold. At present, this is set to use 85% of the bulk bond length as the penalty threshold, where interatomic distances below this are heavily penalised (see https://hiphive.materialsmodeling.org/_modules/hiphive/structure_generation/rattle.html). This choice was motivated by looking at V_Cd and Cd_i as example cases. Note that this is also set to ignore the defect and the distorted neighbours, as these can have bond lengths smaller than this.

It could be a good idea to have a minimum interatomic distance for creating the distortions (as well as for rattling), to avoid certain interstitials where -60/50% gives structures that will explode. Tests are currently underway to see if we can determine a reasonable default minimum interatomic distance to apply with distortions (but can optionally be changed)(note AIRSS typically uses 1.5 Angstrom for this https://airss-docs.github.io/tutorials/examples/)

The penalty for small interatomic distances should help reduce the cases where calculations crash/explode, due to too small initial interatomic distances, and the sporadic cases where the relaxations get stuck in a high-energy metastable minimum.

This is an Issue to track the development of this package, to have it ready for JOSS submission.

Key tasks and the best order of operation (as I see it):

• Docstrings for each function (completed with 982e708)
• Set up GitHub Actions to auto-run the tests when a new push is made.
• Re-modularise. Once the tests and docstrings are in place, we should think about the best module structure for the package, and re-organise the scripts a bit.
• Examples. The current example notebook is a great template, and the test data can also be used for this.
• Tests for the key functions in each module. Currently have this for the distortions.py, input.py, energy_lowering_distortions.py, analysis.py, plotting.py, io.py and cli.py scripts. Most tests now written, just need final scan to check full coverage.
• Feature addition (#3, #4 and any other key features we decide are necessary)
• Documentation (i.e. a docs site). Should be relatively straightforward to generate now that the docstrings are in order, as can be generated from these with just some extra info.

Of course suggestions on this welcome!

In my case, the snb-run command sends jobs to HPS. But the server doesn't execute them (assigned PD or CF state) and after some time an error message appears in the defect directory (see attached).

At the same time, when manually running the task with the sbutch command, the task is calculated.

Could you please tell, where I am using this command incorrectly?

snb-run1_log.txt
slurm-2500867_out.txt
run_relax_sh.txt

Hello,

Using ShakeNBreak for the first time - it's great!

I've just installed on my HPC using conda (conda install --channel conda-forge shakenbreak) in a fresh conda environment.

Running snb-run -h gives the error cannot import name _format_defect_name from doped.utils.plotting'.

I fixed it easy enough by removing the leading underscore from _format_defect_name where it occurs in cli.py, plotting.py and energy_lowering_distortions.py, and can see that your main branch also has these changes.

So I think a new conda release might fix this for other people? Just flagging up. Thanks!

It's been suggested defect-finder sounds too similar to interstitial position locators, so something more focused on the symmetry-breaking / distortion behaviour might be better.

What about ShakeNBreak? As we're shaking the atoms and breaking the symmetry? If you're happy with this @ireaml, then I can rename to this.

greetings! It's Xuecong here and I used the command in 'For a single defect' part of 'Structure Generation' to generate the intrinsic defect of my Cr2O3 calculation.
snb-generate --bulk POSCAR_hex_120.vasp --defect POSCAR
(here .vasp is my 120 symmetry origin bulk and POSCAR has a manully added O interstitial atom)
However I find out that the atom relative far from the defect cell(O1) are shifted even further tham the atom around my defect. Could you help with this problem? Many thanks
following is the POSCAR_hex_120.vasp and POACAR
input.txt

This is an Issue to track the development of a potential feature of a 'localised rattle' function.

I realise the original 'localised rattle' function did work ok, as once the structure was run through the vasp_input.vasp_gam_files() function, it would reorder the POSCAR to be correct (initially appeared as if it was outputting structures in a format that would give the wrong POTCAR/POSCAR ordering).

Either way, the local rattle is not currently recommended or being used by anyone that I know of, so better to reimplement as an improved version, where the rattling displacement amplitude tails off as we move away from the defect site, rather than a hard cutoff radius. To do this, it will likely require a modified version of the hiphive mc_rattle function, with code similar to that here:
https://github.com/materialsproject/pymatgen/blob/834dc9c507f0ceb735887be2db358be37c52a4d7/pymatgen/transformations/site_transformations.py#L584 (similar idea but for spherical distortions).

This would likely avoid the performance decrease witnessed for 'localised rattling' with a hard cutoff radius (which is believed to be due to an initial long-range crystal potential still favouring the high-symmetry defect structure), by still breaking the long-range symmetry, but potentially with the added advantages of being more likely to avoid the sporadic cases of getting stuck in high-energy metastable states, and potentially a small speed boost (but likely marginal, based off previous tests).

Hey @ireaml!
Something I was thinking about again. I wonder if it would be best to choose the rattle stdev based off the bulk bond length. We already get d_min from the bulk bond length, so easy to implement. I know this was briefly discussed before but I can’t really remember why we decided against it 😅
I know one of the issues was that the cases where some defects would consistently return energies higher than unperturbed (due to the supercell being stuck in a higher-energy rattled local minimum) didn’t have a clear correlation as to when it would/wouldn't happen, with it often only occurring for one or two specific defects in specific charge states.

From testing more on a range of compounds here in 🇯🇵, it seems there is definitely a correlation with the bulk bond length, though no clear cutoff or anything. Here I'm plotting the distributions of the bulk bond lengths of oxides I've ran SnB on, for the oxygen vacancies, and coloured orange for any that had a defect species where this occurred.

• I would note here that this was a very rare occurrence (recent updates to setting d_min I think have def helped this), happened in <5% of defects and this is for oxides (which we know are the most likely for this to occur due to the strong ionic bonding), without adjusting the default rattle stdev of 0.25 Å.

From the two datasets (19 and 33 compounds each), you can see there definitely is a trend with this being more likely for the compounds with lower bond lengths:

I know from tests before, we typically found the larger the rattle stdev the better for performance, but I think the effect was typically relatively small, with the key effect of rattling being to break symmetry and disrupt the lattice potential. So do you think it would make sense to set the rattle stdev based on the bulk bond length? Maybe something like 7.5% or 10%?
For comparison we know CdTe and Sb2Se3 were all fine with 0.25 Å, which is around 8.8% and 9.2% each, and I think Sb2S3 didn't have this issue either right? (0.25 Å would be 10% for it). But TiO2 (for which 0.25 Å would be 13%) needed to have it reduced for some cases I think?

If so, I'll rerun each of the cases where this occurred with the rattle stdevs set to 7.5% and 10% of the bulk bond length, and see how they get on! Will also test this choice on the cases I have with GGA that found energy-lowering distortions, and see if using this choice still finds the groundstate

In step 3 there is a bash script to parse all the relaxed energies into files which are used later on in the code. There are a few typos which I have fixed in the below script:

#!/bin/bash

for defect_name in vac_* ; # for each defect
    do cd ${defect_name}/BDM ;
    for i in *; # for each BDM distortion
        do k=$(tail -n1 ${i}/vasp_gam/OUTCAR) ;
        if [[ $k == *"Voluntary"* ]] ;  # check if calc is done
            then echo $i >> ${defect_name}.txt ; # add BDM distortion to txt file
            grep -a sigma ${i}/vasp_gam/OUTCAR | tail -n1 | awk '{print $NF}' >> ${defect_name}.txt ; #and its E
 
        fi ;
        done ;
     cd ../..;
     done

Otherwise works great thank you!

Luisa

This is an Issue to track the feature development of structure I/O compatibility with CP2K and Quantum Espresso.

If possible, this package should work with cp2k, CASTEP, FHI-aims and Quantum Espresso. SIESTA could be a potential addition afterwards (Make note on README that these can be added if users request it?), but less popular and so less of a priority. Should be relatively handy as this just means converting input structures to pymatgen, and then outputs to their format instead of VASP. Likely best implemented as an io.py module?

To do the distortions, we rely on doped reading in the defect type etc., so will have to have some functions that can read an input POSCAR and figure out what defect is present (i.e. the doped parsing tools, which are pymatgen under the hood).

If possible, would be best to not require doped, but have it as an easy way to provide inputs, but alternatively can just give a list of folders in which the defect structures are located, and can have pymatgen read in from these.

The current folder structure generated by ShakeNBreak is defect_name/distortion_type/defect_name_distortion_key/vasp_gam. The structure makes sense, but I think it might be a little clunky / offputting for users.
The vasp_gam part is a legacy effect from the doped vasp_gam_files function, but I think isn't necessary as the user knows these are vasp_gam distortion test runs (from the upper folder structure, and likely they're running the distortion tests in a separate directory to final production run defects anyway). Can be factored out when modifying the VASP I/O parts.
The distortion_type which is either BDM or champion I think can be removed, as again the user should know these are distortion tests from the -XX%_ labels in the distortion_keys etc. For the champion runs, I think these could be made to run in the same defect_name folder as the distortions, just with an obvious different name (like defect_name_champion?)

So suggesting that our default folder structure is changed to: defect_name/defect_name_distortion_key to make it more succinct. What do you think @ireaml?

Hello ,

I have been trying to use the code yet when i am trying to import the library i get the same error

`~/opt/anaconda3/lib/python3.9/site-packages/shakenbreak/distortions.py in
7 import numpy as np
8 from ase.neighborlist import NeighborList
----> 9 from hiphive.structure_generation.rattle import (
10 _probability_mc_rattle,
11 generate_mc_rattled_structures,

~/opt/anaconda3/lib/python3.9/site-packages/hiphive/init.py in
4
5 from .cluster_space import ClusterSpace
...
----> 3 from numba.np.ufunc import _internal
4 from numba.np.ufunc.parallel import ParallelUFuncBuilder, ParallelGUFuncBuilder
5

SystemError: initialization of _internal failed without raising an exception`

I tried changing the numpy and numba version but it doesnt seem to solve the issue,
I was wondering if you could give me some suggestions

The pyproject.toml for shakenbreak 3.3.1 indicates that it relies on doped>=2.3.1, but the most recent PyPI release of doped is 2.3.0.