From a5a8aedad44a2c22bc8524a1e50b6aec24a9901c Mon Sep 17 00:00:00 2001 From: Aliaksandr Yakutovich Date: Tue, 30 May 2023 15:39:04 +0000 Subject: [PATCH 1/4] Update pre-commit config and start applying it. --- .pre-commit-config.yaml | 84 ++- .prospector.yaml | 19 - LICENSE.txt | 2 +- README.md | 2 +- __init__.py | 0 compare.ipynb | 4 +- export_structure.ipynb | 53 +- nanoribbon/editors/__init__.py | 9 +- nanoribbon/editors/replicate.py | 56 +- nanoribbon/viewers/__init__.py | 18 +- nanoribbon/viewers/cdxml2gnr.py | 346 +++++----- nanoribbon/viewers/igor.py | 167 ++--- nanoribbon/viewers/pdos_computed.py | 641 +++++++++---------- nanoribbon/viewers/search.py | 559 ++++++++++------- nanoribbon/viewers/show_computed.py | 901 +++++++++++++++------------ nanoribbon/viewers/smiles2GNR.py | 249 -------- nanoribbon/viewers/smiles2gnr.py | 249 -------- nanoribbon/viewers/structures.py | 685 -------------------- nanoribbon/viewers/utils.py | 82 ++- nanoribbon/workchains/__init__.py | 7 - nanoribbon/workchains/ks_orbitals.py | 255 -------- search.ipynb | 4 +- setup.cfg | 31 +- setup.py | 2 +- start.py | 14 +- submit_ks.ipynb | 124 ---- 26 files changed, 1666 insertions(+), 2897 deletions(-) delete mode 100644 .prospector.yaml delete mode 100644 __init__.py delete mode 100644 nanoribbon/viewers/smiles2GNR.py delete mode 100644 nanoribbon/viewers/smiles2gnr.py delete mode 100644 nanoribbon/viewers/structures.py delete mode 100644 nanoribbon/workchains/__init__.py delete mode 100644 nanoribbon/workchains/ks_orbitals.py delete mode 100644 submit_ks.ipynb diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b68d1ee..602eaaa 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,20 +1,66 @@ -# # Install pre-commit hooks via -# pre-commit install +--- +ci: + autoupdat_schedule: quarterly -- repo: local - hooks: - # yapf = yet another python formatter - - id: yapf - name: yapf - entry: yapf - language: system - types: [python] - args: ["-i"] - - # prospector: collection of linters - - id: prospector - language: system - types: [file, python] - name: prospector - description: "This hook runs Prospector: https://github.com/landscapeio/prospector" - entry: prospector +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.4.0 + hooks: + - id: end-of-file-fixer + - id: trailing-whitespace + - id: check-yaml + - id: check-added-large-files + - repo: https://github.com/pycqa/isort + rev: 5.12.0 + hooks: + - id: isort + args: [--profile, black, --filter-files] + - repo: https://github.com/PyCQA/autoflake + rev: v2.1.1 + hooks: + - id: autoflake + - repo: https://github.com/asottile/pyupgrade + rev: v3.4.0 + hooks: + - id: pyupgrade + args: [--py38-plus] + - repo: https://github.com/psf/black + rev: 23.3.0 + hooks: + - id: black + language_version: python3 # Should be a command that runs python3.6+ + - repo: https://github.com/PyCQA/flake8 + rev: 6.0.0 + hooks: + - id: flake8 + args: [--count, --show-source, --statistics] + additional_dependencies: + - flake8-bugbear + - flake8-builtins + - flake8-comprehensions + - flake8-debugger + - flake8-logging-format + - pep8-naming + - pyflakes + - tryceratops + - repo: https://github.com/pre-commit/mirrors-mypy + rev: v1.3.0 + hooks: + - id: mypy + additional_dependencies: + - types-click-spinner + - types-requests + - types-tabulate + - types-toml + - repo: https://github.com/jumanjihouse/pre-commit-hook-yamlfmt + rev: 0.2.3 + hooks: + - id: yamlfmt + - repo: https://github.com/asottile/setup-cfg-fmt + rev: v2.2.0 + hooks: + - id: setup-cfg-fmt + - repo: https://github.com/kynan/nbstripout + rev: 0.6.1 + hooks: + - id: nbstripout diff --git a/.prospector.yaml b/.prospector.yaml deleted file mode 100644 index 3067378..0000000 --- a/.prospector.yaml +++ /dev/null @@ -1,19 +0,0 @@ -max-line-length: 120 - -ignore-paths: - - doc - - examples - - test - - utils - -pylint: - run: true - -pyflakes: - run: false - -pep8: - run: false - -mccabe: - run: false diff --git a/LICENSE.txt b/LICENSE.txt index a6a7277..dacf60f 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -18,4 +18,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. \ No newline at end of file +SOFTWARE. diff --git a/README.md b/README.md index f089aff..3fe05b7 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ To get an idea of how it works, please [check out the videos](https://www.youtub ## Uploading a new structure File formats compatible with [ASE](https://wiki.fysik.dtu.dk/ase/) can be used. -In the example a .mol file for a structure is generated with ChemDraw. +In the example a .mol file for a structure is generated with ChemDraw. The structure is uploaded, a warning message is provided if similar structures were already computed. diff --git a/__init__.py b/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/compare.ipynb b/compare.ipynb index ab1a80b..caf9175 100644 --- a/compare.ipynb +++ b/compare.ipynb @@ -52,9 +52,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "def make_plots(pk_list):\n", diff --git a/export_structure.ipynb b/export_structure.ipynb index 5e01ca6..12336aa 100644 --- a/export_structure.ipynb +++ b/export_structure.ipynb @@ -2,24 +2,9 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "7834195bf63045f38d46585f64d12e3f", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - "_ColormakerRegistry()" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "%aiida\n", "from aiidalab_widgets_base import viewer\n", @@ -28,24 +13,9 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "b13b606b7c574f279828ee27314fffdf", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - "StructureDataVisualizer(children=(NGLWidget(), HBox(children=(Dropdown(description='File format:', options=('x…" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "uuid = urlparse.parse_qs(urlparse.urlsplit(jupyter_notebook_url).query)['uuid'][0]\n", "display(viewer(load_node(uuid=uuid)))" @@ -53,20 +23,9 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "'676167bc-1bae-4b70-a5cc-d813731f1db2'" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "uuid" ] diff --git a/nanoribbon/editors/__init__.py b/nanoribbon/editors/__init__.py index 4a4bab9..5cd1e9f 100644 --- a/nanoribbon/editors/__init__.py +++ b/nanoribbon/editors/__init__.py @@ -1,8 +1,9 @@ -"""Nanoribbon AiiDA lab viewers.""" -# pylint: disable=unused-import,wrong-import-position +"""Nanoribbon AiiDAlab viewers.""" -from __future__ import absolute_import from aiida import load_profile + load_profile() -from .replicate import NanoribbonReplicateEditor \ No newline at end of file +from .replicate import NanoribbonReplicateEditor + +__all__ = ["NanoribbonReplicateEditor"] diff --git a/nanoribbon/editors/replicate.py b/nanoribbon/editors/replicate.py index dda3884..121ad48 100644 --- a/nanoribbon/editors/replicate.py +++ b/nanoribbon/editors/replicate.py @@ -1,48 +1,34 @@ -import numpy as np -from numpy.linalg import norm -from ase import Atoms -from ase.data import covalent_radii -from ase.neighborlist import NeighborList -import ase.neighborlist -import scipy.stats -from scipy.constants import physical_constants -import itertools -from IPython.display import display, clear_output, HTML -import nglview import ipywidgets as ipw - -from traitlets import HasTraits, Instance, Dict, Unicode, dlink, link, observe -from aiidalab_widgets_base import StructureManagerWidget +from ase import Atoms +from traitlets import Instance class NanoribbonReplicateEditor(ipw.VBox): structure = Instance(Atoms, allow_none=True) - def __init__(self, title=''): + + def __init__(self, title=""): self.title = title self._molecule = None - self.nx_slider = ipw.IntSlider(description="nx", min=1, max=6, continuous_update=False) - + self.nx_slider = ipw.IntSlider( + description="nx", min=1, max=6, continuous_update=False + ) + self.create_bttn = ipw.Button(description="Replicate") self.create_bttn.on_click(self.replicate) - self.info = ipw.HTML('') - super().__init__(children=[ - ipw.HBox([self.nx_slider, self.create_bttn]), - self.info, - ]) - + self.info = ipw.HTML("") + super().__init__( + children=[ + ipw.HBox([self.nx_slider, self.create_bttn]), + self.info, + ] + ) + def replicate(self, _=None): """Create slab and remember the last molecule used.""" - #sa = StructureAnalyzer() - #sa.structure = self.molecule - self.info.value = '' + # sa = StructureAnalyzer() + # sa.structure = self.molecule + self.info.value = "" atoms = self.structure.copy() - nx = self.nx_slider.value - + self.nx_slider.value + self.structure = atoms.repeat((self.nx_slider.value, 1, 1)) - -# @observe('molecule') -# def on_struct_change(self, change=None): -# """Selected molecule from structure.""" -# -# if self.molecule: -# self.nx_slider.value = 1 diff --git a/nanoribbon/viewers/__init__.py b/nanoribbon/viewers/__init__.py index 40bf5c0..4b36fbf 100644 --- a/nanoribbon/viewers/__init__.py +++ b/nanoribbon/viewers/__init__.py @@ -1,12 +1,18 @@ -"""Nanoribbon AiiDA lab viewers.""" -# pylint: disable=unused-import,wrong-import-position +"""Nanoribbon AiiDAlab viewers.""" + -from __future__ import absolute_import from aiida import load_profile + load_profile() +from .cdxml2gnr import CdxmlUpload2GnrWidget +from .pdos_computed import NanoribbonPDOSWidget from .search import NanoribbonSearchWidget from .show_computed import NanoribbonShowWidget -from .pdos_computed import NanoribbonPDOSWidget -#from .smiles2gnr import Smiles2GnrWidget -from .cdxml2gnr import CdxmlUpload2GnrWidget + +__all__ = [ + "CdxmlUpload2GnrWidget", + "NanoribbonPDOSWidget", + "NanoribbonSearchWidget", + "NanoribbonShowWidget", +] diff --git a/nanoribbon/viewers/cdxml2gnr.py b/nanoribbon/viewers/cdxml2gnr.py index 74b68b3..97349af 100644 --- a/nanoribbon/viewers/cdxml2gnr.py +++ b/nanoribbon/viewers/cdxml2gnr.py @@ -1,70 +1,87 @@ # pylint: disable=no-member """Widget to convert SMILES to nanoribbons.""" -import numpy as np -from scipy.stats import mode import re -from IPython.display import clear_output +import ase.neighborlist import ipywidgets as ipw import nglview - -from traitlets import Instance - +import numpy as np from ase import Atoms -from ase.data import covalent_radii,chemical_symbols +from ase.data import chemical_symbols, covalent_radii from ase.neighborlist import NeighborList -import ase.neighborlist - - -from sklearn import manifold, datasets +from IPython.display import clear_output +from scipy.stats import mode from sklearn.decomposition import PCA +from traitlets import Instance + try: import pybel as pb except ModuleNotFoundError: - from openbabel import pybel as pb - + from openbabel import pybel as pb class CdxmlUpload2GnrWidget(ipw.VBox): """Class that allows to upload structures from user's computer.""" + structure = Instance(Atoms, allow_none=True) - def __init__(self, title='CDXML to GNR', description="Upload Structure"): + def __init__(self, title="CDXML to GNR", description="Upload Structure"): try: import openbabel # pylint: disable=unused-import except ImportError: super().__init__( - [ipw.HTML("The SmilesWidget requires the OpenBabel library, " - "but the library was not found.")]) + [ + ipw.HTML( + "The SmilesWidget requires the OpenBabel library, " + "but the library was not found." + ) + ] + ) return - + self.title = title - self.mols=None + self.mols = None self.original_structure = None - self.selection = set() - self.file_upload = ipw.FileUpload(description=description, multiple=False, layout={'width': 'initial'}) + self.selection = set() + self.file_upload = ipw.FileUpload( + description=description, multiple=False, layout={"width": "initial"} + ) supported_formats = ipw.HTML( """ Supported structure formats: ".cdxml" - """) - - self.file_upload.observe(self._on_file_upload, names='value') - - self.create_cell_btn = ipw.Button(description="create GNR", button_style='info', disabled=True) + """ + ) + + self.file_upload.observe(self._on_file_upload, names="value") + + self.create_cell_btn = ipw.Button( + description="create GNR", button_style="info", disabled=True + ) self.create_cell_btn.on_click(self._on_cell_button_pressed) self.viewer = nglview.NGLWidget() - self.viewer.stage.set_parameters(mouse_preset='pymol') - self.viewer.observe(self._on_picked, names='picked') - self.allmols = ipw.Dropdown(options=[None],description='Select mol',value=None, disabled=True) - self.allmols.observe(self._on_sketch_selected,names='value') + self.viewer.stage.set_parameters(mouse_preset="pymol") + self.viewer.observe(self._on_picked, names="picked") + self.allmols = ipw.Dropdown( + options=[None], description="Select mol", value=None, disabled=True + ) + self.allmols.observe(self._on_sketch_selected, names="value") self.select_two = ipw.HTML("") self.picked_out = ipw.Output() - self.cell_button_out = ipw.Output() - super().__init__(children=[self.file_upload, supported_formats,self.allmols,self.select_two, - self.viewer, self.picked_out, self.create_cell_btn, self.cell_button_out]) - + self.cell_button_out = ipw.Output() + super().__init__( + children=[ + self.file_upload, + supported_formats, + self.allmols, + self.select_two, + self.viewer, + self.picked_out, + self.create_cell_btn, + self.cell_button_out, + ] + ) @staticmethod def guess_scaling_factor(atoms): @@ -94,45 +111,51 @@ def guess_scaling_factor(atoms): # Scale box to match equilibrium carbon-carbon bond distance. cc_eq = 1.4313333333 - return cc_eq / avg_bond - - @staticmethod + return cc_eq / avg_bond + + @staticmethod def scale(atoms, s): """Scale atomic positions by the `factor`.""" c_x, c_y, c_z = atoms.cell atoms.set_cell((s * c_x, s * c_y, c_z), scale_atoms=True) atoms.center() return atoms - - @staticmethod - def pybel2ase(mol): + + @staticmethod + def pybel2ase(mol): """converts pybel molecule into ase Atoms""" - asemol = Atoms() - species=[chemical_symbols[atm.atomicnum] for atm in mol.atoms] - pos=np.asarray([atm.coords for atm in mol.atoms]) + Atoms() + species = [chemical_symbols[atm.atomicnum] for atm in mol.atoms] + pos = np.asarray([atm.coords for atm in mol.atoms]) pca = PCA(n_components=3) - posnew=pca.fit_transform(pos) + posnew = pca.fit_transform(pos) atoms = Atoms(species, positions=posnew) - sys_size = np.ptp(atoms.positions,axis=0) - atoms.rotate(-90, 'z') #cdxml are rotated - atoms.pbc=True + sys_size = np.ptp(atoms.positions, axis=0) + atoms.rotate(-90, "z") # cdxml are rotated + atoms.pbc = True atoms.cell = sys_size + 10 atoms.center() return atoms - @staticmethod + @staticmethod def construct_cell(atoms, id1, id2): """Construct periodic cell based on two selected equivalent atoms.""" - pos = [[atoms[id2].x, atoms[id2].y], [atoms[id1].x, atoms[id1].y], [atoms[id2].x + int(1), atoms[id1].y]] + pos = [ + [atoms[id2].x, atoms[id2].y], + [atoms[id1].x, atoms[id1].y], + [atoms[id2].x + int(1), atoms[id1].y], + ] vec = [np.array(pos[0]) - np.array(pos[1]), np.array(pos[2]) - np.array(pos[1])] c_x = np.linalg.norm(vec[0]) angle = np.math.atan2(np.linalg.det([vec[0], vec[1]]), np.dot(vec[0], vec[1])) if np.abs(angle) > 0.01: - atoms.rotate_euler(center=atoms[id1].position, phi=-angle, theta=0.0, psi=0.0) + atoms.rotate_euler( + center=atoms[id1].position, phi=-angle, theta=0.0, psi=0.0 + ) c_y = 15.0 + np.amax(atoms.positions[:, 1]) - np.amin(atoms.positions[:, 1]) c_z = 15.0 + np.amax(atoms.positions[:, 2]) - np.amin(atoms.positions[:, 2]) @@ -141,85 +164,91 @@ def construct_cell(atoms, id1, id2): atoms.pbc = (True, True, True) atoms.center() atoms.wrap(eps=0.001) - -# #### REMOVE EXTRA ATOMS -# some_too_close=True -# while some_too_close and len(atoms)>1: -# #ite+=1 -# some_too_close=False -# cov_radii = [covalent_radii[a.number] for a in atoms] -# nl = NeighborList(cov_radii, bothways = True, self_interaction = False) -# nl.update(atoms) -# for a in atoms: -# if some_too_close: -# break -# -# indices, offsets = nl.get_neighbors(a.index) -# #print(ite,a,indices) -# for i, offset in zip(indices, offsets): -# dist = np.linalg.norm(a.position -(atoms.positions[i] + np.dot(offset, atoms.get_cell()))) -# if dist < 0.4 : -# some_too_close=True -# #remove H atoms before others -# if atoms[i].symbol != 'H' and a.symbol == 'H': -# del atoms[a.index] -# #print('deleted ',a.index) -# else: -# del atoms[i] -# #print('deleted ',i) -# break -# -# elif a.symbol == "C" : #this is executed only if dist > 0.4 -# -# #this part removes H from overcoordinated C -# #if atoms[i].symbol=='H' and [ida.symbol for ida in atoms[indices]].count('C') == 3: -# -# #this part removes H from all C and N -# if atoms[i].symbol=='H' or a.symbol=='N' : -# some_too_close=True -# del atoms[i] -# #print('deleted ',i) -# break -# -# -# #### ADD Hydrogens -# cov_radii = [covalent_radii[a.number] for a in atoms] -# nl = NeighborList(cov_radii, bothways = True, self_interaction = False) -# nl.update(atoms) -# -# need_a_H = [] -# for a in atoms: -# nlist=nl.get_neighbors(a.index)[0] -# if len(nlist)<3: -# if a.symbol=='C' or a.symbol=='N': -# need_a_H.append(a.index) -# -# #print("Added missing Hydrogen atoms: ", need_a_H) -# -# dCH=1.1 -# for a in need_a_H: -# vec = np.zeros(3) -# indices, offsets = nl.get_neighbors(atoms[a].index) -# for i, offset in zip(indices, offsets): -# vec += -atoms[a].position +(atoms.positions[i] + np.dot(offset, atoms.get_cell())) -# vec = -vec/np.linalg.norm(vec)*dCH -# vec += atoms[a].position -# htoadd = ase.Atom('H',vec) -# atoms.append(htoadd) -# -# return atoms - + + # #### REMOVE EXTRA ATOMS + # some_too_close=True + # while some_too_close and len(atoms)>1: + # #ite+=1 + # some_too_close=False + # cov_radii = [covalent_radii[a.number] for a in atoms] + # nl = NeighborList(cov_radii, bothways = True, self_interaction = False) + # nl.update(atoms) + # for a in atoms: + # if some_too_close: + # break + # + # indices, offsets = nl.get_neighbors(a.index) + # #print(ite,a,indices) + # for i, offset in zip(indices, offsets): + # dist = np.linalg.norm(a.position -(atoms.positions[i] + np.dot(offset, atoms.get_cell()))) + # if dist < 0.4 : + # some_too_close=True + # #remove H atoms before others + # if atoms[i].symbol != 'H' and a.symbol == 'H': + # del atoms[a.index] + # #print('deleted ',a.index) + # else: + # del atoms[i] + # #print('deleted ',i) + # break + # + # elif a.symbol == "C" : #this is executed only if dist > 0.4 + # + # #this part removes H from overcoordinated C + # #if atoms[i].symbol=='H' and [ida.symbol for ida in atoms[indices]].count('C') == 3: + # + # #this part removes H from all C and N + # if atoms[i].symbol=='H' or a.symbol=='N' : + # some_too_close=True + # del atoms[i] + # #print('deleted ',i) + # break + # + # + # #### ADD Hydrogens + # cov_radii = [covalent_radii[a.number] for a in atoms] + # nl = NeighborList(cov_radii, bothways = True, self_interaction = False) + # nl.update(atoms) + # + # need_a_H = [] + # for a in atoms: + # nlist=nl.get_neighbors(a.index)[0] + # if len(nlist)<3: + # if a.symbol=='C' or a.symbol=='N': + # need_a_H.append(a.index) + # + # #print("Added missing Hydrogen atoms: ", need_a_H) + # + # dCH=1.1 + # for a in need_a_H: + # vec = np.zeros(3) + # indices, offsets = nl.get_neighbors(atoms[a].index) + # for i, offset in zip(indices, offsets): + # vec += -atoms[a].position +(atoms.positions[i] + np.dot(offset, atoms.get_cell())) + # vec = -vec/np.linalg.norm(vec)*dCH + # vec += atoms[a].position + # htoadd = ase.Atom('H',vec) + # atoms.append(htoadd) + # + # return atoms # Remove redundant atoms. ORIGINAL tobedel = [] - n_l = NeighborList([covalent_radii[a.number] for a in atoms], bothways=False, self_interaction=False) + n_l = NeighborList( + [covalent_radii[a.number] for a in atoms], + bothways=False, + self_interaction=False, + ) n_l.update(atoms) for atm in atoms: indices, offsets = n_l.get_neighbors(atm.index) for i, offset in zip(indices, offsets): - dist = np.linalg.norm(atm.position - (atoms.positions[i] + np.dot(offset, atoms.get_cell()))) + dist = np.linalg.norm( + atm.position + - (atoms.positions[i] + np.dot(offset, atoms.get_cell())) + ) if dist < 0.4: tobedel.append(atoms[i].index) @@ -228,13 +257,17 @@ def construct_cell(atoms, id1, id2): # Find unit cell and apply it. ## Add Hydrogens. - n_l = NeighborList([covalent_radii[a.number] for a in atoms], bothways=True, self_interaction=False) + n_l = NeighborList( + [covalent_radii[a.number] for a in atoms], + bothways=True, + self_interaction=False, + ) n_l.update(atoms) need_hydrogen = [] for atm in atoms: if len(n_l.get_neighbors(atm.index)[0]) < 3: - if atm.symbol == 'C' or atm.symbol=='N' : + if atm.symbol == "C" or atm.symbol == "N": need_hydrogen.append(atm.index) print("Added missing Hydrogen atoms: ", need_hydrogen) @@ -243,16 +276,18 @@ def construct_cell(atoms, id1, id2): vec = np.zeros(3) indices, offsets = n_l.get_neighbors(atoms[atm].index) for i, offset in zip(indices, offsets): - vec += -atoms[atm].position + (atoms.positions[i] + np.dot(offset, atoms.get_cell())) + vec += -atoms[atm].position + ( + atoms.positions[i] + np.dot(offset, atoms.get_cell()) + ) vec = -vec / np.linalg.norm(vec) * 1.1 + atoms[atm].position - atoms.append(ase.Atom('H', vec)) + atoms.append(ase.Atom("H", vec)) return atoms def _on_picked(self, _=None): """When an attom is picked.""" - if 'atom1' not in self.viewer.picked.keys(): + if "atom1" not in self.viewer.picked.keys(): return # did not click on atom self.create_cell_btn.disabled = True @@ -263,7 +298,7 @@ def _on_picked(self, _=None): self.viewer.component_0.remove_ball_and_stick() self.viewer.add_ball_and_stick() - idx = self.viewer.picked['atom1']['index'] + idx = self.viewer.picked["atom1"]["index"] # Toggle. if idx in self.selection: @@ -274,43 +309,54 @@ def _on_picked(self, _=None): if len(self.selection) == 2: self.create_cell_btn.disabled = False - #if(selection): + # if(selection): sel_str = ",".join([str(i) for i in sorted(self.selection)]) print("Selected atoms: " + sel_str) - self.viewer.add_representation('ball+stick', selection="@" + sel_str, color='red', aspectRatio=3.0) - self.viewer.picked = {} # reset, otherwise immidiately selecting same atom again won't create change event + self.viewer.add_representation( + "ball+stick", selection="@" + sel_str, color="red", aspectRatio=3.0 + ) + self.viewer.picked = ( + {} + ) # reset, otherwise immidiately selecting same atom again won't create change event - def _on_file_upload(self, change=None): """When file upload button is pressed.""" - self.mols=None + self.mols = None self.create_cell_btn.disabled = True listmols = [] molid = 0 - for fname, item in change['new'].items(): - frmt = fname.split('.')[-1] - if frmt == 'cdxml': - cdxml_file_string = self.file_upload.value[fname]['content'].decode('ascii') - self.mols=re.findall('') + m = pb.readstring("cdxml", "") self.mols[molid] = m - listmols.append((str(molid)+': '+m.formula, molid)) ## m MUST BE a pb object!!! + listmols.append( + (str(molid) + ": " + m.formula, molid) + ) ## m MUST BE a pb object!!! molid += 1 - self.allmols.options = listmols + self.allmols.options = listmols self.allmols.value = 0 - self.allmols.disabled=False - break - - def _on_sketch_selected(self,change=None): - self.structure = None #needed to empty view in second viewer + self.allmols.disabled = False + break + + def _on_sketch_selected(self, change=None): + self.structure = None # needed to empty view in second viewer if self.mols is None or self.allmols.value is None: return self.create_cell_btn.disabled = True atoms = self.pybel2ase(self.mols[self.allmols.value]) factor = self.guess_scaling_factor(atoms) struct = self.scale(atoms, factor) - self.select_two.value = '

Select two equivalent atoms that define the basis vector

' + self.select_two.value = ( + "

Select two equivalent atoms that define the basis vector

" + ) self.original_structure = struct.copy() if hasattr(self.viewer, "component_0"): self.viewer.component_0.remove_ball_and_stick() @@ -324,7 +370,7 @@ def _on_sketch_selected(self,change=None): self.viewer.add_component(nglview.ASEStructure(struct)) # adds ball+stick self.viewer.center() self.viewer.handle_resize() - self.file_upload.value.clear() + self.file_upload.value.clear() def _on_cell_button_pressed(self, _=None): """Generate GNR button pressed.""" @@ -340,15 +386,17 @@ def _on_cell_button_pressed(self, _=None): incoming_struct = self.original_structure.copy() self.structure = self.construct_cell(self.original_structure, id1, id2) self.original_structure = incoming_struct.copy() - + if hasattr(self.viewer, "component_0"): self.viewer.component_0.remove_ball_and_stick() cid = self.viewer.component_0.id - self.viewer.remove_component(cid) + self.viewer.remove_component(cid) # Empty selection. self.selection = set() # Add new component. - self.viewer.add_component(nglview.ASEStructure(self.original_structure)) # adds ball+stick + self.viewer.add_component( + nglview.ASEStructure(self.original_structure) + ) # adds ball+stick self.viewer.center() - self.viewer.handle_resize() + self.viewer.handle_resize() diff --git a/nanoribbon/viewers/igor.py b/nanoribbon/viewers/igor.py index 9bd258e..d50f623 100644 --- a/nanoribbon/viewers/igor.py +++ b/nanoribbon/viewers/igor.py @@ -5,9 +5,11 @@ """ import re + import numpy as np -class Axis(object): + +class Axis: """Represents an axis of an IGOR wave""" def __init__(self, symbol, min, delta, unit, wavename=None): @@ -22,18 +24,24 @@ def __str__(self): Note: SetScale/P expects minimum value and step-size """ delta = 0 if self.delta is None else self.delta - s = "X SetScale/P {symb} {min},{delta}, \"{unit}\", {name};\n"\ - .format(symb=self.symbol, min=self.min, delta=delta,\ - unit=self.unit, name=self.wavename) + s = 'X SetScale/P {symb} {min},{delta}, "{unit}", {name};\n'.format( + symb=self.symbol, + min=self.min, + delta=delta, + unit=self.unit, + name=self.wavename, + ) return s def read(self, string): """Read axis from string - Format: + Format: X SetScale/P x 0,2.01342281879195e-11,"m", data_00381_Up; SetScale d 0,0,"V", data_00381_Up """ - match = re.search("SetScale/?P? (.) ([+-\.\de]+),([+-\.\de]+),\"(\w+)\",\s*(\w+)", string) + match = re.search( + 'SetScale/?P? (.) ([+-\\.\\de]+),([+-\\.\\de]+),"(\\w+)",\\s*(\\w+)', string + ) self.symbol = match.group(1) self.min = float(match.group(2)) self.delta = float(match.group(3)) @@ -41,8 +49,7 @@ def read(self, string): self.wavename = match.group(5) - -class Wave(object): +class Wave: """A class template for IGOR waves of generic dimension""" def __init__(self, data, axes, name=None): @@ -58,10 +65,10 @@ def __str__(self): dimstring = "(" for i in range(len(self.data.shape)): - dimstring += "{}, ".format(self.data.shape[i]) - dimstring = dimstring[:-2] + ")" + dimstring += f"{self.data.shape[i]}, " + dimstring = dimstring[:-2] + ")" - s += "WAVES/N={} {}\n".format(dimstring, self.name) + s += f"WAVES/N={dimstring} {self.name}\n" s += "BEGIN\n" s += self.print_data() s += "END\n" @@ -71,36 +78,36 @@ def __str__(self): def read(self, fname): """Read IGOR wave - + Should work for any dimension. Tested so far only for 2d wave. """ - f=open(fname, 'r') - content=f.read() + f = open(fname) + content = f.read() f.close() lines = content.split("\r") line = lines.pop(0) if not line == "IGOR": - raise IOError("Files does not begin with 'IGOR'") + raise OSError("Files does not begin with 'IGOR'") line = lines.pop(0) - while not re.match("WAVES",line): + while not re.match("WAVES", line): line = lines.pop(0) - match = re.search("WAVES/N=\(([\d,]+)\)\s+(.+)",line) - grid = match.group(1).split(',') + match = re.search(r"WAVES/N=\(([\d,]+)\)\s+(.+)", line) + grid = match.group(1).split(",") grid = np.array(grid, dtype=int) self.name = match.group(2) line = lines.pop(0) if not line == "BEGIN": - raise IOError("Missing 'BEGIN' statement of data block") + raise OSError("Missing 'BEGIN' statement of data block") # read data datastring = "" line = lines.pop(0) - while not re.match("END",line): + while not re.match("END", line): datastring += line line = lines.pop(0) data = np.array(datastring.split(), dtype=float) @@ -111,13 +118,13 @@ def read(self, fname): matches = re.findall("SetScale.+?(?:;|$)", line) self.axes = [] for match in matches: - ax = Axis(None,None,None,None) + ax = Axis(None, None, None, None) ax.read(match) self.axes.append(ax) # the rest is discarded... - #line = lines.pop(0) - #print(line) + # line = lines.pop(0) + # print(line) @property def extent(self): @@ -127,118 +134,130 @@ def extent(self): for i in range(len(grid)): ax = self.axes[i] extent.append(ax.min) - extent.append(ax.min+ax.delta*grid[i]) + extent.append(ax.min + ax.delta * grid[i]) return np.array(extent) - def print_data(self): """Determines how to print the data block. - + To be implemented by subclasses.""" def write(self, fname): - f=open(fname, 'w') + f = open(fname, "w") f.write(str(self)) f.close() class Wave1d(Wave): """1d Igor wave""" - + default_parameters = dict( - xmin = 0.0, - xdelta = None, - xlabel = 'x', - ylabel = 'y', + xmin=0.0, + xdelta=None, + xlabel="x", + ylabel="y", ) def __init__(self, data=None, axes=None, name="1d", **kwargs): """Initialize 1d IGOR wave""" - super(Wave1d, self).__init__(data, axes, name) + super().__init__(data, axes, name) self.parameters = self.default_parameters for key, value in kwargs.items(): if key in self.parameters: self.parameters[key] = value else: - raise KeyError("Unknown parameter {}".format(key)) + raise KeyError(f"Unknown parameter {key}") if axes is None: - p=self.parameters - x = Axis(symbol='x', min=p['xmin'], delta=p['xdelta'], unit=p['xlabel'], - wavename=self.name) + p = self.parameters + x = Axis( + symbol="x", + min=p["xmin"], + delta=p["xdelta"], + unit=p["xlabel"], + wavename=self.name, + ) self.axes = [x] def print_data(self): s = "" for line in self.data: - s += "{:12.6e}\n".format(float(line)) + s += f"{float(line):12.6e}\n" return s - class Wave2d(Wave): """2d Igor wave""" default_parameters = dict( - xmin = 0.0, - xdelta = None, - xmax = None, - xlabel = 'x', - ymin = 0.0, - ydelta = None, - ymax = None, - ylabel = 'y', + xmin=0.0, + xdelta=None, + xmax=None, + xlabel="x", + ymin=0.0, + ydelta=None, + ymax=None, + ylabel="y", ) - + def __init__(self, data=None, axes=None, name=None, **kwargs): """Initialize 2d Igor wave Parameters ---------- - - * data - * name - * xmin, xdelta, xlabel - * ymin, ydelta, ylabel + + * data + * name + * xmin, xdelta, xlabel + * ymin, ydelta, ylabel """ - super(Wave2d, self).__init__(data, axes=axes, name=name) + super().__init__(data, axes=axes, name=name) self.parameters = self.default_parameters for key, value in kwargs.items(): if key in self.parameters: self.parameters[key] = value else: - raise KeyError("Unknown parameter {}".format(key)) + raise KeyError(f"Unknown parameter {key}") if axes is None: - p=self.parameters + p = self.parameters nx, ny = self.data.shape - if p['xmax'] is None: - p['xmax'] = p['xdelta'] * nx - elif p['xdelta'] is None: - p['xdelta'] = p['xmax'] / nx - - if p['ymax'] is None: - p['ymax'] = p['ydelta'] * ny - elif p['ydelta'] is None: - p['ydelta'] = p['ymax'] / ny - - x = Axis(symbol='x', min=p['xmin'], delta=p['xdelta'], - unit=p['xlabel'], wavename=self.name) - y = Axis(symbol='y', min=p['ymin'], delta=p['ydelta'], - unit=p['ylabel'], wavename=self.name) - self.axes = [x,y] - + if p["xmax"] is None: + p["xmax"] = p["xdelta"] * nx + elif p["xdelta"] is None: + p["xdelta"] = p["xmax"] / nx + + if p["ymax"] is None: + p["ymax"] = p["ydelta"] * ny + elif p["ydelta"] is None: + p["ydelta"] = p["ymax"] / ny + + x = Axis( + symbol="x", + min=p["xmin"], + delta=p["xdelta"], + unit=p["xlabel"], + wavename=self.name, + ) + y = Axis( + symbol="y", + min=p["ymin"], + delta=p["ydelta"], + unit=p["ylabel"], + wavename=self.name, + ) + self.axes = [x, y] def print_data(self): """Determines how to print the data block""" s = "" for line in self.data: for x in line: - s += "{:12.6e} ".format(x) + s += f"{x:12.6e} " s += "\n" - return s \ No newline at end of file + return s diff --git a/nanoribbon/viewers/pdos_computed.py b/nanoribbon/viewers/pdos_computed.py index 479e86e..4d8fb79 100644 --- a/nanoribbon/viewers/pdos_computed.py +++ b/nanoribbon/viewers/pdos_computed.py @@ -1,114 +1,98 @@ -import gzip -import re - -import nglview -import gzip -import bqplot as bq -import numpy as np import io +import re from base64 import b64encode from xml.etree import ElementTree -import matplotlib.pyplot as plt -from collections import OrderedDict -import scipy.constants as const -import ase -from ase.data import covalent_radii, atomic_numbers -from ase.data.colors import cpk_colors -from ase.neighborlist import NeighborList -import ase.io -import ase.io.cube -from aiida.orm import CalcJobNode, load_node, QueryBuilder, WorkChainNode - -import scipy.ndimage - -import matplotlib.pyplot as plt -import matplotlib -from matplotlib.collections import LineCollection -from matplotlib.ticker import FormatStrFormatter import ipywidgets as ipw +import matplotlib +import matplotlib.pyplot as plt +import nglview +import numpy as np +import scipy.ndimage +from aiida.orm import CalcJobNode, QueryBuilder, WorkChainNode +from IPython.core.display import HTML from IPython.display import clear_output -from IPython.core.display import HTML, Javascript +from matplotlib.ticker import FormatStrFormatter on_band_click_global = None - - - def get_calc_by_label(workcalc, label): calcs = get_calcs_by_label(workcalc, label) assert len(calcs) == 1 return calcs[0] - + + def get_calcs_by_label(workcalc, label): qb = QueryBuilder() - qb.append(WorkChainNode, filters={'uuid':workcalc.uuid}) - qb.append(CalcJobNode, with_incoming=WorkChainNode, filters={'label':label}) - calcs = [ c[0] for c in qb.all() ] - for calc in calcs: - assert(calc.is_finished_ok == True) - return calcs - -def w0gauss(x,n): + qb.append(WorkChainNode, filters={"uuid": workcalc.uuid}) + qb.append(CalcJobNode, with_incoming=WorkChainNode, filters={"label": label}) + calcs = [c[0] for c in qb.all()] + for calc in calcs: + assert calc.is_finished_ok == True + return calcs + + +def w0gauss(x, n): arg = np.minimum(200.0, x**2) - w0gauss = np.exp ( - arg) / np.sqrt(np.pi) - if n==0 : + w0gauss = np.exp(-arg) / np.sqrt(np.pi) + if n == 0: return w0gauss hd = 0.0 - hp = np.exp( - arg) + hp = np.exp(-arg) ni = 0 a = 1.0 / np.sqrt(np.pi) - for i in range(1, n+1): + for i in range(1, n + 1): hd = 2.0 * x * hp - 2.0 * ni * hd ni = ni + 1 - a = - a / (i * 4.0) - hp = 2.0 * x * hd-2.0 * ni * hp + a = -a / (i * 4.0) + hp = 2.0 * x * hd - 2.0 * ni * hp ni = ni + 1 w0gauss = w0gauss + a * hp return w0gauss + def var_width_lines(x, y, lw, aspect): - nx = len(x) edge_up = [np.zeros(nx), np.zeros(nx)] edge_down = [np.zeros(nx), np.zeros(nx)] - + for i_x in range(nx): if i_x == 0: - dx = x[i_x] - x[i_x+1] - dy = y[i_x] - y[i_x+1] + dx = x[i_x] - x[i_x + 1] + dy = y[i_x] - y[i_x + 1] elif i_x == nx - 1: - dx = x[i_x-1] - x[i_x] - dy = y[i_x-1] - y[i_x] + dx = x[i_x - 1] - x[i_x] + dy = y[i_x - 1] - y[i_x] else: - dx = x[i_x-1] - x[i_x+1] - dy = y[i_x-1] - y[i_x+1] - - line_dir = np.array((dx, dy)) + dx = x[i_x - 1] - x[i_x + 1] + dy = y[i_x - 1] - y[i_x + 1] + + line_dir = np.array((dx, dy)) # Convert line vector to "figure coordinates" line_dir[0] /= aspect - + perp_dir = np.array((line_dir[1], -line_dir[0])) - shift_vec = perp_dir/np.sqrt(perp_dir[0]**2 + perp_dir[1]**2)*lw[i_x] - + shift_vec = perp_dir / np.sqrt(perp_dir[0] ** 2 + perp_dir[1] ** 2) * lw[i_x] + # convert shift_vec back to "data coordinates" shift_vec[0] *= aspect - + edge_up[0][i_x] = x[i_x] + shift_vec[0] edge_up[1][i_x] = y[i_x] + shift_vec[1] edge_down[0][i_x] = x[i_x] - shift_vec[0] edge_down[1][i_x] = y[i_x] - shift_vec[1] - + return edge_up, edge_down + class NanoribbonPDOSWidget(ipw.VBox): def __init__(self, workcalc, **kwargs): self._workcalc = workcalc - self.vacuum_level = workcalc.get_extra('vacuum_level') - self.homo = workcalc.get_extra('homo') - self.lumo = workcalc.get_extra('lumo') + self.vacuum_level = workcalc.get_extra("vacuum_level") + self.homo = workcalc.get_extra("homo") + self.lumo = workcalc.get_extra("lumo") pdos_calc = get_calc_by_label(workcalc, "export_pdos") bands_calc = get_calc_by_label(workcalc, "bands") self.structure = bands_calc.inputs.structure @@ -117,209 +101,274 @@ def __init__(self, workcalc, **kwargs): self.selected_spin = None self.bands = bands_calc.outputs.output_band.get_bands() if self.bands.ndim == 2: - self.bands = self.bands[None,:,:] - - -#1 + self.bands = self.bands[None, :, :] + + # 1 atomic_proj_xml = pdos_calc.outputs.retrieved.open("atomic_proj.xml").name root = ElementTree.parse(atomic_proj_xml).getroot() - if 'NUMBER_OF_BANDS' in root.find('HEADER').attrib: + if "NUMBER_OF_BANDS" in root.find("HEADER").attrib: self.parse_atomic_proj_xml(root) else: self.parse_old_atomic_proj_xml(root) - -#2 - output_log = pdos_calc.outputs.retrieved.open('aiida.out').name + + # 2 + output_log = pdos_calc.outputs.retrieved.open("aiida.out").name # parse mapping atomic functions -> atoms # example: state # 2: atom 1 (C ), wfc 2 (l=1 m= 1) content = open(output_log).read() - m = re.findall("\n\s+state #\s*(\d+): atom\s*(\d+) ", content, re.DOTALL) - self.atmwfc2atom = dict([(int(i), int(j)) for i,j in m]) + m = re.findall("\n\\s+state #\\s*(\\d+): atom\\s*(\\d+) ", content, re.DOTALL) + self.atmwfc2atom = {int(i): int(j) for i, j in m} assert len(self.atmwfc2atom) == self.natwfcs - assert len(set(self.atmwfc2atom.values())) == self.natoms - - -#3 + assert len(set(self.atmwfc2atom.values())) == self.natoms + + # 3 self.kpts = np.linspace(0.0, 0.5, self.nkpoints) - #atmwfcsues, projections = correct_band_crossings(kpts, eigvalues, projections) + # atmwfcsues, projections = correct_band_crossings(kpts, eigvalues, projections) self.bands = np.swapaxes(self.eigvalues, 1, 2) + self.vacuum_level + ####BOH -####BOH - - style = {"description_width":"200px"} + style = {"description_width": "200px"} layout = ipw.Layout(width="600px") - self.sigma_slider = ipw.FloatSlider(description="Broadening [eV]", min=0.01, max=0.5, value=0.1, step=0.01, - continuous_update=False, layout=layout, style=style) - #sigma_slider.observe(on_change, names='value') - self.ngauss_slider = ipw.IntSlider(description="Methfessel-Paxton order", min=0, max=3, value=0, - continuous_update=False, layout=layout, style=style) - #ngauss_slider.observe(on_change, names='value') - - self.colorpicker = ipw.ColorPicker(concise=True, description='PDOS color', value='orange', style=style) - - center = (self.homo + self.lumo)/2.0 - self.emin_box = ipw.FloatText(description="Emin (eV)", value=np.round(center-3.0, 1), step=0.1, style=style) - self.emax_box = ipw.FloatText(description="Emax (eV)", value=np.round(center+3.0, 1), step=0.1, style=style) - self.band_proj_box = ipw.FloatText(description="Max band width (eV)", value=0.1, step=0.01, style=style) - - + self.sigma_slider = ipw.FloatSlider( + description="Broadening [eV]", + min=0.01, + max=0.5, + value=0.1, + step=0.01, + continuous_update=False, + layout=layout, + style=style, + ) + # sigma_slider.observe(on_change, names='value') + self.ngauss_slider = ipw.IntSlider( + description="Methfessel-Paxton order", + min=0, + max=3, + value=0, + continuous_update=False, + layout=layout, + style=style, + ) + # ngauss_slider.observe(on_change, names='value') + + self.colorpicker = ipw.ColorPicker( + concise=True, description="PDOS color", value="orange", style=style + ) + + center = (self.homo + self.lumo) / 2.0 + self.emin_box = ipw.FloatText( + description="Emin (eV)", + value=np.round(center - 3.0, 1), + step=0.1, + style=style, + ) + self.emax_box = ipw.FloatText( + description="Emax (eV)", + value=np.round(center + 3.0, 1), + step=0.1, + style=style, + ) + self.band_proj_box = ipw.FloatText( + description="Max band width (eV)", value=0.1, step=0.01, style=style + ) + + # on_change(None) - #on_change(None) - self.plot_button = ipw.Button(description="Plot") self.plot_button.on_click(self.on_plot_click) - self.selected_atoms = set() + self.selected_atoms = set() self.viewer = nglview.NGLWidget() - self.viewer.add_component(nglview.ASEStructure(self.ase_struct)) # adds ball+stick + self.viewer.add_component( + nglview.ASEStructure(self.ase_struct) + ) # adds ball+stick self.viewer.add_unitcell() self.viewer.center() - - self.viewer.stage.set_parameters(mouse_preset='pymol') - #self.viewer.stage.set_parameters(mouse_preset='coot') - self.viewer.observe(self.on_picked, names='picked') + self.viewer.stage.set_parameters(mouse_preset="pymol") + # self.viewer.stage.set_parameters(mouse_preset='coot') + + self.viewer.observe(self.on_picked, names="picked") self.plot_out = ipw.Output() - boxes=[self.sigma_slider, self.ngauss_slider, self.viewer, - self.emin_box, self.emax_box, self.band_proj_box, - self.colorpicker, self.plot_button, self.plot_out] - - super(NanoribbonPDOSWidget, self).__init__(boxes, **kwargs) - - - + boxes = [ + self.sigma_slider, + self.ngauss_slider, + self.viewer, + self.emin_box, + self.emax_box, + self.band_proj_box, + self.colorpicker, + self.plot_button, + self.plot_out, + ] + + super().__init__(boxes, **kwargs) + def parse_atomic_proj_xml(self, root): - - header = root.find('HEADER') - - self.nbands = int(header.attrib['NUMBER_OF_BANDS']) - self.nkpoints = int(header.attrib['NUMBER_OF_K-POINTS']) - self.nspins = int(header.attrib['NUMBER_OF_SPIN_COMPONENTS']) - self.natwfcs = int(header.attrib['NUMBER_OF_ATOMIC_WFC']) - - + header = root.find("HEADER") + + self.nbands = int(header.attrib["NUMBER_OF_BANDS"]) + self.nkpoints = int(header.attrib["NUMBER_OF_K-POINTS"]) + self.nspins = int(header.attrib["NUMBER_OF_SPIN_COMPONENTS"]) + self.natwfcs = int(header.attrib["NUMBER_OF_ATOMIC_WFC"]) + self.eigvalues = np.zeros((self.nspins, self.nbands, self.nkpoints)) self.kpoint_weights = [] - - self.projections = np.zeros((self.nspins, self.nbands, self.nkpoints, self.natwfcs)) - - eigenstates = root.find('EIGENSTATES') - + + self.projections = np.zeros( + (self.nspins, self.nbands, self.nkpoints, self.natwfcs) + ) + + eigenstates = root.find("EIGENSTATES") + i_spin = 0 i_kpt = -1 - + for child in eigenstates: - - if child.tag == 'K-POINT': + if child.tag == "K-POINT": i_kpt += 1 if i_kpt >= self.nkpoints: i_spin = 1 i_kpt = 0 if i_spin == 0: - self.kpoint_weights.append(float(child.attrib['Weight'])) - - if child.tag == 'E': - - arr = np.fromstring(child.text, sep='\n') - self.eigvalues[i_spin, :, i_kpt] = arr * 13.60569806589 - self.vacuum_level # convert Ry to eV - - if child.tag == 'PROJS': + self.kpoint_weights.append(float(child.attrib["Weight"])) + + if child.tag == "E": + arr = np.fromstring(child.text, sep="\n") + self.eigvalues[i_spin, :, i_kpt] = ( + arr * 13.60569806589 - self.vacuum_level + ) # convert Ry to eV + + if child.tag == "PROJS": for i, c in enumerate(child): - - if c.tag == 'ATOMIC_WFC': - + if c.tag == "ATOMIC_WFC": arr = np.fromstring(c.text.replace(",", "\n"), sep="\n") - arr2 = arr.reshape(self.nbands, 2) # group real and imaginary part together - arr3 = np.sum(np.square(arr2), axis=1) # calculate square of abs value + arr2 = arr.reshape( + self.nbands, 2 + ) # group real and imaginary part together + arr3 = np.sum( + np.square(arr2), axis=1 + ) # calculate square of abs value self.projections[i_spin, :, i_kpt, i] = arr3 - - self.kpoint_weights = np.array(self.kpoint_weights) + self.kpoint_weights = np.array(self.kpoint_weights) def parse_old_atomic_proj_xml(self, root): - - self.nbands = int(root.find('HEADER/NUMBER_OF_BANDS').text) - self.nkpoints = int(root.find('HEADER/NUMBER_OF_K-POINTS').text) - self.nspins = int(root.find('HEADER/NUMBER_OF_SPIN_COMPONENTS').text) - self.natwfcs = int(root.find('HEADER/NUMBER_OF_ATOMIC_WFC').text) + self.nbands = int(root.find("HEADER/NUMBER_OF_BANDS").text) + self.nkpoints = int(root.find("HEADER/NUMBER_OF_K-POINTS").text) + self.nspins = int(root.find("HEADER/NUMBER_OF_SPIN_COMPONENTS").text) + self.natwfcs = int(root.find("HEADER/NUMBER_OF_ATOMIC_WFC").text) - self.kpoint_weights = np.fromstring(root.find('WEIGHT_OF_K-POINTS').text, sep=' ') + self.kpoint_weights = np.fromstring( + root.find("WEIGHT_OF_K-POINTS").text, sep=" " + ) self.eigvalues = np.zeros((self.nspins, self.nbands, self.nkpoints)) for i in range(self.nspins): for k in range(self.nkpoints): - eigtag = 'EIG.%s'%(i+1) if self.nspins > 1 else 'EIG' - arr = np.fromstring(root.find('EIGENVALUES/K-POINT.%d/%s'%(k+1, eigtag)).text, sep='\n') - self.eigvalues[i, :, k] = arr * 13.60569806589 - self.vacuum_level # convert Ry to eV - - self.projections = np.zeros((self.nspins, self.nbands, self.nkpoints, self.natwfcs)) + eigtag = "EIG.%s" % (i + 1) if self.nspins > 1 else "EIG" + arr = np.fromstring( + root.find("EIGENVALUES/K-POINT.%d/%s" % (k + 1, eigtag)).text, + sep="\n", + ) + self.eigvalues[i, :, k] = ( + arr * 13.60569806589 - self.vacuum_level + ) # convert Ry to eV + + self.projections = np.zeros( + (self.nspins, self.nbands, self.nkpoints, self.natwfcs) + ) for i in range(self.nspins): for k in range(self.nkpoints): for l in range(self.natwfcs): - spintag = 'SPIN.%d/'%(i+1) if self.nspins > 1 else "" - raw = root.find('PROJECTIONS/K-POINT.%d/%sATMWFC.%d'%(k+1, spintag, l+1)).text + spintag = "SPIN.%d/" % (i + 1) if self.nspins > 1 else "" + raw = root.find( + "PROJECTIONS/K-POINT.%d/%sATMWFC.%d" % (k + 1, spintag, l + 1) + ).text arr = np.fromstring(raw.replace(",", "\n"), sep="\n") - arr2 = arr.reshape(self.nbands, 2) # group real and imaginary part together - arr3 = np.sum(np.square(arr2), axis=1) # calculate square of abs value + arr2 = arr.reshape( + self.nbands, 2 + ) # group real and imaginary part together + arr3 = np.sum( + np.square(arr2), axis=1 + ) # calculate square of abs value self.projections[i, :, k, l] = arr3 - -#4 - def calc_pdos(self,sigma, ngauss, Emin, Emax, atmwfcs=None): + # 4 + + def calc_pdos(self, sigma, ngauss, Emin, Emax, atmwfcs=None): DeltaE = 0.01 - x = np.arange(Emin,Emax,DeltaE) + x = np.arange(Emin, Emax, DeltaE) # calculate histogram for all spins, bands, and kpoints in parallel - xx = np.tile(x[:, None, None, None], (1, self.nspins, self.nbands, self.nkpoints)) + xx = np.tile( + x[:, None, None, None], (1, self.nspins, self.nbands, self.nkpoints) + ) arg = (xx - self.eigvalues) / sigma delta = w0gauss(arg, n=ngauss) / sigma if atmwfcs: - p = np.sum(self.projections[:,:,:,atmwfcs], axis=3) # sum over selected atmwfcs + p = np.sum( + self.projections[:, :, :, atmwfcs], axis=3 + ) # sum over selected atmwfcs else: - p = np.sum(self.projections, axis=3) # sum over all atmwfcs + p = np.sum(self.projections, axis=3) # sum over all atmwfcs c = delta * p * self.kpoint_weights - y = np.sum(c, axis=(2,3)) # sum over bands and kpoints + y = np.sum(c, axis=(2, 3)) # sum over bands and kpoints return x, y -#5 + # 5 def igor_pdos(self): - center = (self.homo + self.lumo)/2.0 - Emin, Emax = center-3.0, center+3.0 + center = (self.homo + self.lumo) / 2.0 + Emin, Emax = center - 3.0, center + 3.0 if self.selected_atoms: - atmwfcs = [k-1 for k, v in self.atmwfc2atom.items() if v-1 in self.selected_atoms] + atmwfcs = [ + k - 1 + for k, v in self.atmwfc2atom.items() + if v - 1 in self.selected_atoms + ] else: atmwfcs = None - pdos = self.calc_pdos(ngauss=self.ngauss_slider.value, sigma=self.sigma_slider.value, Emin=Emin, Emax=Emax, atmwfcs=atmwfcs) + pdos = self.calc_pdos( + ngauss=self.ngauss_slider.value, + sigma=self.sigma_slider.value, + Emin=Emin, + Emax=Emax, + atmwfcs=atmwfcs, + ) e = pdos[0] p = pdos[1].transpose()[0] tempio = io.StringIO() with tempio as f: - f.write(u'IGOR\rWAVES\te1\tp1\rBEGIN\r') + f.write("IGOR\rWAVES\te1\tp1\rBEGIN\r") for x, y in zip(e, p): - f.write(u'\t{:.8f}\t{:.8f}\r'.format(x, y)) - f.write(u'END\r') - f.write(u'X SetScale/P x 0,1,"", e1; SetScale y 0,0,"", e1\rX SetScale/P x 0,1,"", p1; SetScale y 0,0,"", p1\r') + f.write(f"\t{x:.8f}\t{y:.8f}\r") + f.write("END\r") + f.write( + 'X SetScale/P x 0,1,"", e1; SetScale y 0,0,"", e1\rX SetScale/P x 0,1,"", p1; SetScale y 0,0,"", p1\r' + ) return f.getvalue() def mk_igor_link(self): igorvalue = self.igor_pdos() igorfile = b64encode(igorvalue.encode()).decode() - filename = self.ase_struct.get_chemical_formula() + "_pk%d.itx" % self.structure.pk + filename = ( + self.ase_struct.get_chemical_formula() + "_pk%d.itx" % self.structure.pk + ) - html = ' 0): - return int(num + .5) + if num > 0: + return int(num + 0.5) else: - return int(num - .5) + return int(num - 0.5) if fermi_energy and number_electrons: - raise ValueError("Specify either the number of electrons or the " - "Fermi energy, but not both") + raise ValueError( + "Specify either the number of electrons or the " + "Fermi energy, but not both" + ) - assert bandsdata.units == 'eV' + assert bandsdata.units == "eV" stored_bands = bandsdata.get_bands() if len(stored_bands.shape) == 3: @@ -284,15 +364,16 @@ def nint(num): # analysis on occupations: if fermi_energy is None: - num_kpoints = len(bands) if number_electrons is None: try: _, stored_occupations = bandsdata.get_bands(also_occupations=True) except KeyError: - raise KeyError("Cannot determine metallicity if I don't have " - "either fermi energy, or occupations") + raise KeyError( + "Cannot determine metallicity if I don't have " + "either fermi energy, or occupations" + ) # put the occupations in the same order of bands, also in case of multiple bands if len(stored_occupations.shape) == 3: @@ -300,7 +381,9 @@ def nint(num): # spin up and spin down array # put all spins on one band per kpoint - occupations = np.concatenate([_ for _ in stored_occupations], axis=1) + occupations = np.concatenate( + [_ for _ in stored_occupations], axis=1 + ) else: occupations = stored_occupations @@ -309,22 +392,42 @@ def nint(num): # sort the bands by energy, and reorder the occupations accordingly # since after joining the two spins, I might have unsorted stuff - bands, occupations = [np.array(y) for y in zip(*[zip(*j) for j in - [sorted(zip(i[0].tolist(), i[1].tolist()), - key=lambda x: x[0]) - for i in zip(bands, occupations)]])] - number_electrons = int(round(sum([sum(i) for i in occupations]) / num_kpoints)) - - homo_indexes = [np.where(np.array([nint(_) for _ in x]) > 0)[0][-1] for x in occupations] - if len(set(homo_indexes)) > 1: # there must be intersections of valence and conduction bands + bands, occupations = ( + np.array(y) + for y in zip( + *[ + zip(*j) + for j in [ + sorted( + zip(i[0].tolist(), i[1].tolist()), + key=lambda x: x[0], + ) + for i in zip(bands, occupations) + ] + ] + ) + ) + number_electrons = int( + round(sum([sum(i) for i in occupations]) / num_kpoints) + ) + + homo_indexes = [ + np.where(np.array([nint(_) for _ in x]) > 0)[0][-1] + for x in occupations + ] + if ( + len(set(homo_indexes)) > 1 + ): # there must be intersections of valence and conduction bands return False, None, None, None else: homo = [_[0][_[1]] for _ in zip(bands, homo_indexes)] try: lumo = [_[0][_[1] + 1] for _ in zip(bands, homo_indexes)] except IndexError: - raise ValueError("To understand if it is a metal or insulator, " - "need more bands than n_band=number_electrons") + raise ValueError( + "To understand if it is a metal or insulator, " + "need more bands than n_band=number_electrons" + ) else: bands = np.sort(bands) @@ -334,13 +437,19 @@ def nint(num): # calculation, 2 otherwise) number_electrons_per_band = 4 - len(stored_bands.shape) # 1 or 2 # gather the energies of the homo band, for every kpoint - homo = [i[number_electrons / number_electrons_per_band - 1] for i in bands] # take the nth level + homo = [ + i[number_electrons / number_electrons_per_band - 1] for i in bands + ] # take the nth level try: # gather the energies of the lumo band, for every kpoint - lumo = [i[number_electrons / number_electrons_per_band] for i in bands] # take the n+1th level + lumo = [ + i[number_electrons / number_electrons_per_band] for i in bands + ] # take the n+1th level except IndexError: - raise ValueError("To understand if it is a metal or insulator, " - "need more bands than n_band=number_electrons") + raise ValueError( + "To understand if it is a metal or insulator, " + "need more bands than n_band=number_electrons" + ) if number_electrons % 2 == 1 and len(stored_bands.shape) == 2: # if #electrons is odd and we have a non spin polarized calculation @@ -349,9 +458,9 @@ def nint(num): # if the nth band crosses the (n+1)th, it is an insulator gap = min(lumo) - max(homo) - if gap == 0.: - return False, 0., None, None - elif gap < 0.: + if gap == 0.0: + return False, 0.0, None, None + elif gap < 0.0: return False, gap, None, None else: return True, gap, max(homo), min(lumo) @@ -367,21 +476,26 @@ def nint(num): max_mins = [(max(i), min(i)) for i in levels] if fermi_energy > bands.max(): - raise ValueError("The Fermi energy is above all band energies, " - "don't know what to do") + raise ValueError( + "The Fermi energy is above all band energies, " + "don't know what to do" + ) if fermi_energy < bands.min(): - raise ValueError("The Fermi energy is below all band energies, " - "don't know what to do.") + raise ValueError( + "The Fermi energy is below all band energies, " + "don't know what to do." + ) # one band is crossed by the fermi energy if any(i[1] < fermi_energy and fermi_energy < i[0] for i in max_mins): - return False, 0., None, None + return False, 0.0, None, None # case of semimetals, fermi energy at the crossing of two bands # this will only work if the dirac point is computed! - elif (any(i[0] == fermi_energy for i in max_mins) and - any(i[1] == fermi_energy for i in max_mins)): - return False, 0., None, None + elif any(i[0] == fermi_energy for i in max_mins) and any( + i[1] == fermi_energy for i in max_mins + ): + return False, 0.0, None, None # insulating case else: # take the max of the band maxima below the fermi energy @@ -390,9 +504,10 @@ def nint(num): lumo = min([i[1] for i in max_mins if i[1] > fermi_energy]) gap = lumo - homo - if gap <= 0.: - raise Exception("Something wrong has been implemented. " - "Revise the code!") + if gap <= 0.0: + raise Exception( + "Something wrong has been implemented. " "Revise the code!" + ) return True, gap, homo, lumo def search(self, do_all=False): @@ -405,47 +520,47 @@ def search(self, do_all=False): self.results.value = "searching..." # html table header - html = '' + html = "" html += '
' html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' - html += '' + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" + html += "" # query AiiDA database filters = {} - filters['label'] = 'NanoribbonWorkChain' - filters['extras.preprocess_version'] = self.PREPROCESS_VERSION - filters['extras.preprocess_successful'] = True - filters['extras.obsolete'] = False + filters["label"] = "NanoribbonWorkChain" + filters["extras.preprocess_version"] = self.PREPROCESS_VERSION + filters["extras.preprocess_successful"] = True + filters["extras.obsolete"] = False pk_list = self.inp_pks.value.strip().split() if pk_list: - filters['id'] = {'in': pk_list} + filters["id"] = {"in": pk_list} formula_list = self.inp_formula.value.strip().split() if self.inp_formula.value: - filters['extras.formula'] = {'in': formula_list} + filters["extras.formula"] = {"in": formula_list} if len(self.text_description.value) > 1: - filters['description'] = {'like': '%{}%'.format(text_description.value)} + filters["description"] = {"like": f"%{text_description.value}%"} def add_range_filter(bounds, label): - filters['extras.'+label] = {'and':[{'>=':bounds[0]}, {'<':bounds[1]}]} + filters["extras." + label] = {"and": [{">=": bounds[0]}, {"<": bounds[1]}]} add_range_filter(self.inp_gap.value, "gap") add_range_filter(self.inp_homo.value, "homo") @@ -454,55 +569,61 @@ def add_range_filter(bounds, label): add_range_filter(self.inp_tmagn.value, "total_magnetization_per_angstr") add_range_filter(self.inp_amagn.value, "absolute_magnetization_per_angstr") - qb = QueryBuilder() + qb = QueryBuilder() qb.append(WorkChainNode, filters=filters) - qb.order_by({WorkChainNode:{'ctime':'desc'}}) + qb.order_by({WorkChainNode: {"ctime": "desc"}}) for i, node_tuple in enumerate(qb.iterall()): node = node_tuple[0] - thumbnail = node.get_extra('thumbnail') - description = node.get_extra('structure_description') - opt_structure_uuid = node.get_extra('opt_structure_uuid') - opt_structure_pk = node.get_extra('opt_structure_pk') + thumbnail = node.get_extra("thumbnail") + node.get_extra("structure_description") + opt_structure_uuid = node.get_extra("opt_structure_uuid") + opt_structure_pk = node.get_extra("opt_structure_pk") # append table row - html += '' - html += ''.format(node.id) - html += ''.format(node.pk, node.pk) - html += '' % node.ctime.strftime("%Y-%m-%d %H:%M") - html += '' % node.get_extra('formula') - html += '' % node.description - html += '' % node.get_extra('homo') - html += '' % node.get_extra('lumo') - html += '' % node.get_extra('gap') - html += '' % node.get_extra('fermi_energy') - html += '' % node.get_extra('energy') - html += '' % node.get_extra('cellx') - html += '' % node.get_extra('total_magnetization_per_angstr') - html += '' % node.get_extra('absolute_magnetization_per_angstr') - html += '' % (node.inputs['tot_charge'] if 'tot_charge' in node.inputs else 0.0) - html += ''.format(thumbnail, opt_structure_pk) - html += ''.format(node.id) - html += '' - - html += '
PKCreatedFormulaCalcNameHOMO (eV)LUMO (eV)GAP (eV)Fermi (eV)Energy (eV)Cell x (Å)Total Mag./ÅAbs Mag./ÅChargeStructure
PKCreatedFormulaCalcNameHOMO (eV)LUMO (eV)GAP (eV)Fermi (eV)Energy (eV)Cell x (Å)Total Mag./ÅAbs Mag./ÅChargeStructure
{}%s%s%s%4.2f%4.2f%4.2f%4.2f%9.3f%5.2f%4.2f%4.2f%3.1f'.format(opt_structure_uuid) - html += 'Show
'.format(node.id) - html += 'PDOS
' - html += 'Found {} matching entries.
'.format(qb.count()) + html += "" + html += f'' + html += f'{node.pk}' + html += "%s" % node.ctime.strftime("%Y-%m-%d %H:%M") + html += "%s" % node.get_extra("formula") + html += "%s" % node.description + html += "%4.2f" % node.get_extra("homo") + html += "%4.2f" % node.get_extra("lumo") + html += "%4.2f" % node.get_extra("gap") + html += "%4.2f" % node.get_extra("fermi_energy") + html += "%9.3f" % node.get_extra("energy") + html += "%5.2f" % node.get_extra("cellx") + html += "%4.2f" % node.get_extra("total_magnetization_per_angstr") + html += "%4.2f" % node.get_extra( + "absolute_magnetization_per_angstr" + ) + html += "%3.1f" % ( + node.inputs["tot_charge"] if "tot_charge" in node.inputs else 0.0 + ) + html += f'' + html += f'' + html += ( + f'Show
' + ) + html += f'PDOS' + html += "" + + html += "" + html += f"Found {qb.count()} matching entries.
" html += '' - html += '
' + html += "" self.results.value = html def update_filter_limits(self): """Query the database to initalize the extremal values of the filters.""" filters = {} - filters['label'] = 'NanoribbonWorkChain' - filters['extras.preprocess_version'] = self.PREPROCESS_VERSION - filters['extras.preprocess_successful'] = True - filters['extras.obsolete'] = False + filters["label"] = "NanoribbonWorkChain" + filters["extras.preprocess_version"] = self.PREPROCESS_VERSION + filters["extras.preprocess_successful"] = True + filters["extras.obsolete"] = False - qb = QueryBuilder() + qb = QueryBuilder() qb.append(WorkChainNode, filters=filters) def compare_and_set(val, slider): @@ -515,9 +636,13 @@ def compare_and_set(val, slider): for m in qb.all(): calc = m[0] - compare_and_set(calc.get_extra('gap'), self.inp_gap) - compare_and_set(calc.get_extra('homo'), self.inp_homo) - compare_and_set(calc.get_extra('lumo'), self.inp_lumo) - compare_and_set(calc.get_extra('fermi_energy'), self.inp_efermi) - compare_and_set(calc.get_extra('total_magnetization_per_angstr'), self.inp_tmagn) - compare_and_set(calc.get_extra('absolute_magnetization_per_angstr'), self.inp_amagn) + compare_and_set(calc.get_extra("gap"), self.inp_gap) + compare_and_set(calc.get_extra("homo"), self.inp_homo) + compare_and_set(calc.get_extra("lumo"), self.inp_lumo) + compare_and_set(calc.get_extra("fermi_energy"), self.inp_efermi) + compare_and_set( + calc.get_extra("total_magnetization_per_angstr"), self.inp_tmagn + ) + compare_and_set( + calc.get_extra("absolute_magnetization_per_angstr"), self.inp_amagn + ) diff --git a/nanoribbon/viewers/show_computed.py b/nanoribbon/viewers/show_computed.py index 4edc0aa..2ee0d47 100644 --- a/nanoribbon/viewers/show_computed.py +++ b/nanoribbon/viewers/show_computed.py @@ -1,35 +1,40 @@ -# -*- coding: utf-8 -*- """Viewers to display the results of the Nanoribbon work chain.""" +import copy + # Base imports. import gzip -import re import io import tempfile from base64 import b64encode from collections import OrderedDict -import nglview -import ipywidgets as ipw + +import ase +import ase.io.cube import bqplot as bq -import numpy as np -import matplotlib.pyplot as plt +import ipywidgets as ipw import matplotlib -from IPython.display import display, clear_output +import matplotlib.pyplot as plt +import nglview +import numpy as np import scipy.constants as const -import ase -import ase.io.cube -from traitlets import dlink, observe, Instance, Int - -import copy # AiiDA imports. from aiida.common import exceptions from aiida.plugins import DataFactory +from IPython.display import clear_output, display +from traitlets import Instance, Int, dlink, observe -# Local imports. -from .utils import plot_struct_2d, get_calc_by_label, get_calcs_by_label, from_cube_to_arraydata from .igor import Wave2d +# Local imports. +from .utils import ( + from_cube_to_arraydata, + get_calc_by_label, + get_calcs_by_label, + plot_struct_2d, +) + # AiiDA data objects. ArrayData = DataFactory("array") # pylint: disable=invalid-name BandsData = DataFactory("array.bands") # pylint: disable=invalid-name @@ -48,62 +53,75 @@ class BandsViewerWidget(ipw.VBox): """Widget to view AiiDA BandsData object.""" + bands = Instance(BandsData, allow_none=True) structure = Instance(StructureData, allow_none=True) selected_band = Int(allow_none=True) selected_kpoint = Int(allow_none=True) selected_spin = Int(allow_none=True) selected_3D = Int(allow_none=True) - + def __init__(self, **kwargs): - self.bands = kwargs['bands'] - self.structure = kwargs['structure'] - self.vacuum_level = kwargs['vacuum_level'] + self.bands = kwargs["bands"] + self.structure = kwargs["structure"] + self.vacuum_level = kwargs["vacuum_level"] self.bands_array = self.bands.get_bands() self.band_plots = [] - self.homo = kwargs['homo'] - self.lumo = kwargs['lumo'] - - self.num_export_bands = kwargs['num_export_bands'] + self.homo = kwargs["homo"] + self.lumo = kwargs["lumo"] + + self.num_export_bands = kwargs["num_export_bands"] # Always make the array 3-dimensional. - nptk_ks = 12 #this is hardcoded in bands_lowres for the moment + nptk_ks = 12 # this is hardcoded in bands_lowres for the moment if self.bands_array.ndim == 2: self.bands_array = self.bands_array[None, :, :] self.eff_mass_parabolas = [] layout = ipw.Layout(padding="5px", margin="0px") - layout = ipw.Layout(padding="5px", margin="0px", width='auto') + layout = ipw.Layout(padding="5px", margin="0px", width="auto") # Slider to control how many points of the band to use for parabolic fit. - self.efm_fit_slider = ipw.IntSlider(description="Eff. mass fit", - min=3, - max=15, - step=2, - continuous_update=False, - layout=layout) - band_selector = ipw.IntSlider(description="Band", - value=int(kwargs['nelectrons'] / 2) , - min=max(1,int(kwargs['nelectrons'] / 2) - self.num_export_bands//2), - max=int(kwargs['nelectrons'] / 2) + self.num_export_bands//2, - step=1, - continuous_update=False, - readout_format='d', - layout=layout) - kpoint_slider = ipw.IntSlider(description="k-point", - min=1, - max=nptk_ks, - readout_format='d', - continuous_update=False, - layout=layout) - spin_selector = ipw.RadioButtons(options=[('up', 0), ('down', 1)], - description='Select spin', - disabled=False) - view_3D = ipw.RadioButtons(options=[('no', 0), ('yes', 1)], - description='plot3D', - disabled=False) - - boxes = [self.efm_fit_slider, band_selector, kpoint_slider, spin_selector,view_3D] + self.efm_fit_slider = ipw.IntSlider( + description="Eff. mass fit", + min=3, + max=15, + step=2, + continuous_update=False, + layout=layout, + ) + band_selector = ipw.IntSlider( + description="Band", + value=int(kwargs["nelectrons"] / 2), + min=max(1, int(kwargs["nelectrons"] / 2) - self.num_export_bands // 2), + max=int(kwargs["nelectrons"] / 2) + self.num_export_bands // 2, + step=1, + continuous_update=False, + readout_format="d", + layout=layout, + ) + kpoint_slider = ipw.IntSlider( + description="k-point", + min=1, + max=nptk_ks, + readout_format="d", + continuous_update=False, + layout=layout, + ) + spin_selector = ipw.RadioButtons( + options=[("up", 0), ("down", 1)], description="Select spin", disabled=False + ) + view_3D = ipw.RadioButtons( + options=[("no", 0), ("yes", 1)], description="plot3D", disabled=False + ) + + boxes = [ + self.efm_fit_slider, + band_selector, + kpoint_slider, + spin_selector, + view_3D, + ] plots = [] for ispin in range(self.bands_array.shape[0]): @@ -113,10 +131,10 @@ def __init__(self, **kwargs): self.eff_mass_parabolas.append(eff_mass_parabola) boxes.append(ipw.HBox(plots)) - dlink((kpoint_slider, 'value'), (self, 'selected_kpoint')) - dlink((band_selector, 'value'), (self, 'selected_band')) - dlink((spin_selector, 'value'), (self, 'selected_spin')) - dlink((view_3D, 'value'), (self, 'selected_3D')) + dlink((kpoint_slider, "value"), (self, "selected_kpoint")) + dlink((band_selector, "value"), (self, "selected_band")) + dlink((spin_selector, "value"), (self, "selected_spin")) + dlink((view_3D, "value"), (self, "selected_3D")) # Display the orbital map also initially. self.on_band_change(_=None) @@ -138,70 +156,81 @@ def plot_bands(self, ispin): x_data = np.linspace(0.0, x_max, nkpoints) y_datas = self.bands_array[ispin, :, :].transpose() - self.vacuum_level - lines = bq.Lines(x=x_data, - y=y_datas, - color=np.zeros(nbands), - animate=True, - stroke_width=4.0, - scales={ - 'x': x_sc, - 'y': y_sc, - 'color': bq.ColorScale(colors=['gray', 'red'], min=0.0, max=1.0) - }) - - homo_line = bq.Lines(x=[0, x_max], - y=[center, center], - line_style='dashed', - colors=['red'], - scales={ - 'x': x_sc, - 'y': y_sc - }) + lines = bq.Lines( + x=x_data, + y=y_datas, + color=np.zeros(nbands), + animate=True, + stroke_width=4.0, + scales={ + "x": x_sc, + "y": y_sc, + "color": bq.ColorScale(colors=["gray", "red"], min=0.0, max=1.0), + }, + ) + + homo_line = bq.Lines( + x=[0, x_max], + y=[center, center], + line_style="dashed", + colors=["red"], + scales={"x": x_sc, "y": y_sc}, + ) # Initialize the parabola as a random line and set visible to false # Later, when it is correctly set, show it. - eff_mass_parabola = bq.Lines(x=[0, 0], - y=[0, 0], - visible=False, - stroke_width=1.0, - line_style='solid', - colors=['blue'], - scales={ - 'x': x_sc, - 'y': y_sc - }) - ax_x = bq.Axis(label=u'kA^-1', scale=x_sc, grid_lines='solid', tick_format='.3f', tick_values=[0, x_max]) - ax_y = bq.Axis(label='eV', scale=y_sc, orientation='vertical', grid_lines='solid') - - fig = bq.Figure(axes=[ax_x, ax_y], - marks=[lines, homo_line, eff_mass_parabola], - title='Spin {}'.format(ispin), - layout=ipw.Layout(height="800px", width="200px"), - fig_margin={ - "left": 45, - "top": 60, - "bottom": 60, - "right": 40 - }, - min_aspect_ratio=0.25, - max_aspect_ratio=0.25) + eff_mass_parabola = bq.Lines( + x=[0, 0], + y=[0, 0], + visible=False, + stroke_width=1.0, + line_style="solid", + colors=["blue"], + scales={"x": x_sc, "y": y_sc}, + ) + ax_x = bq.Axis( + label="kA^-1", + scale=x_sc, + grid_lines="solid", + tick_format=".3f", + tick_values=[0, x_max], + ) + ax_y = bq.Axis( + label="eV", scale=y_sc, orientation="vertical", grid_lines="solid" + ) + + fig = bq.Figure( + axes=[ax_x, ax_y], + marks=[lines, homo_line, eff_mass_parabola], + title=f"Spin {ispin}", + layout=ipw.Layout(height="800px", width="200px"), + fig_margin={"left": 45, "top": 60, "bottom": 60, "right": 40}, + min_aspect_ratio=0.25, + max_aspect_ratio=0.25, + ) save_btn = ipw.Button(description="Download png") - save_btn.on_click(lambda b: fig.save_png()) # save_png() does not work with unicode labels + save_btn.on_click( + lambda b: fig.save_png() + ) # save_png() does not work with unicode labels - box = ipw.VBox([fig, save_btn, self.mk_igor_link(ispin)], - layout=ipw.Layout(align_items="center", padding="5px", margin="0px")) + box = ipw.VBox( + [fig, save_btn, self.mk_igor_link(ispin)], + layout=ipw.Layout(align_items="center", padding="5px", margin="0px"), + ) return box, lines, eff_mass_parabola def mk_igor_link(self, ispin): """Create a downloadable link.""" igorvalue = self.igor_bands(ispin) igorfile = b64encode(igorvalue.encode()).decode() - filename = self.structure.get_ase().get_chemical_formula() + "_bands_spin{}_pk{}.itx".format( - ispin, self.structure.id) + filename = ( + self.structure.get_ase().get_chemical_formula() + + f"_bands_spin{ispin}_pk{self.structure.id}.itx" + ) - html = ' 0: vmin = 0.0 + + if vmax < 0: + vmax = 0.0 + if vmin > 0: + vmin = 0.0 x_2 = self._current_structure.cell[0][0] * 2.0 y_2 = self._current_structure.cell[1][1] # Set labels and limits. - axplt.set_xlabel(u'Å') - axplt.set_ylabel(u'Å') + axplt.set_xlabel("Å") + axplt.set_ylabel("Å") axplt.set_xlim(0, x_2) axplt.set_ylim(0, y_2) - + if self.center0: amax = np.max(np.abs([vmax, vmin])) vmin = -amax vmax = amax - - plot = axplt.imshow(np.tile(flipped_data, (1, 2)), extent=[0, x_2, 0, y_2], - cmap=self.cmap, vmin=vmin, vmax=vmax) - + + plot = axplt.imshow( + np.tile(flipped_data, (1, 2)), + extent=[0, x_2, 0, y_2], + cmap=self.cmap, + vmin=vmin, + vmax=vmax, + ) + if self.show_cbar: - fig.colorbar(plot, - label=self.units, - ticks=[vmin, vmax], - format='%.e', - orientation='horizontal', - shrink=0.3) + fig.colorbar( + plot, + label=self.units, + ticks=[vmin, vmax], + format="%.e", + orientation="horizontal", + shrink=0.3, + ) # To show structure uniformly with transparency, add fully opaque structure and # then an additional transparent data layer on top plot_struct_2d(axplt, self._current_structure, 1.0) - axplt.imshow(np.tile(flipped_data, (1, 2)), extent=[0, x_2, 0, y_2], - cmap=self.cmap, vmin=vmin, vmax=vmax, alpha=(1-self.opacity_slider.value), zorder=10) - + axplt.imshow( + np.tile(flipped_data, (1, 2)), + extent=[0, x_2, 0, y_2], + cmap=self.cmap, + vmin=vmin, + vmax=vmax, + alpha=(1 - self.opacity_slider.value), + zorder=10, + ) + plt.show() - + self.make_export_links(data) - + def make_export_links(self, data): html = "    Download: " html += self.make_export_link_txt(data) html += " " html += self.make_export_link_igor(data) - + self.dl_link.value = html - + def make_export_link_txt(self, data): - header = "xlim=(%.2f, %.2f), ylim=(%.2f, %.2f)" % (0.0, self._current_structure.cell[0][0], - 0.0, self._current_structure.cell[1][1]) + header = "xlim=({:.2f}, {:.2f}), ylim=({:.2f}, {:.2f})".format( + 0.0, + self._current_structure.cell[0][0], + 0.0, + self._current_structure.cell[1][1], + ) tempio = io.BytesIO() np.savetxt(tempio, data, header=header, fmt="%.4e") - + enc_file = b64encode(tempio.getvalue()).decode() - + if self.export_label is not None: plot_name = self.export_label else: plot_name = "export" - - filename = "{}_h{:.3f}.txt".format(plot_name, self.selected_height) - html = 'Spin density")]) + spin_density_vbox.children += tuple( + [ipw.HTML(value="

Spin density

")] + ) spin_density_vbox.children += tuple([self.spin_view_3D]) spin_density_vbox.children += tuple([self.output_s]) - + self.on_spin_view_mode_change() - + # --------------------------------------- # STS mapping widget self.sts_heading = ipw.HTML(value="

LDOS mappings

") - self.sts_energy_text = ipw.FloatText(value=round(self._workcalc.get_extra('homo'), 2), description='energy [eV]:') - self.sts_fwhm_text = ipw.FloatText(value=0.1, description='fwhm [eV]:') - self.sts_mapping_viewer = CubeArrayData2dViewerWidget(cmap='gist_heat', center0=False, - show_cbar=False, export_label="pk%d_ldos" % workcalc.pk) + self.sts_energy_text = ipw.FloatText( + value=round(self._workcalc.get_extra("homo"), 2), description="energy [eV]:" + ) + self.sts_fwhm_text = ipw.FloatText(value=0.1, description="fwhm [eV]:") + self.sts_mapping_viewer = CubeArrayData2dViewerWidget( + cmap="gist_heat", + center0=False, + show_cbar=False, + export_label="pk%d_ldos" % workcalc.pk, + ) self.sts_mapping_viewer_wrapper = ipw.VBox([]) - self.sts_btn = ipw.Button(description='view LDOS') + self.sts_btn = ipw.Button(description="view LDOS") self.sts_btn.on_click(self.on_sts_btn_click) - - self.sts_viewer_box = ipw.VBox([ - self.sts_heading, - ipw.HBox([self.sts_energy_text, self.sts_fwhm_text, self.sts_btn]), - self.sts_mapping_viewer_wrapper - ]) + + self.sts_viewer_box = ipw.VBox( + [ + self.sts_heading, + ipw.HBox([self.sts_energy_text, self.sts_fwhm_text, self.sts_btn]), + self.sts_mapping_viewer_wrapper, + ] + ) # --------------------------------------- - + self.output = ipw.Output() self.on_kpoint_change() - - super().__init__([ - self.info, - ipw.HBox([self.bands_viewer, self.output]), spin_density_vbox, self.sts_viewer_box - ], **kwargs) - + super().__init__( + [ + self.info, + ipw.HBox([self.bands_viewer, self.output]), + spin_density_vbox, + self.sts_viewer_box, + ], + **kwargs, + ) + def on_spin_view_mode_change(self, _=None): - def _set_viewer_cube_data(viewer): try: viewer.arraydata = self.spindensity_calc.outputs.output_data except exceptions.NotExistent: print("Compatibility mode for old spin cube file") - with gzip.open(self.spindensity_calc.outputs.retrieved.open("_spin.cube.gz").name) as fpointer: + with gzip.open( + self.spindensity_calc.outputs.retrieved.open("_spin.cube.gz").name + ) as fpointer: arrayd = from_cube_to_arraydata(fpointer.read()) viewer.arraydata = arrayd - + with self.output_s: clear_output() if self.spin_view_3D.value == 0: @@ -720,65 +842,76 @@ def _set_viewer_cube_data(viewer): else: display(self.spinden_viewer_3d) _set_viewer_cube_data(self.spinden_viewer_3d) - - + def on_view_3D_change(self, _=None): """Plot 3D orbitals in case of selection.""" if self.bands_viewer.selected_3D: - self.twod_3D=[self.orbital_viewer_2d, self.orbital_viewer_3d] - else: - self.twod_3D=[self.orbital_viewer_2d] + self.twod_3D = [self.orbital_viewer_2d, self.orbital_viewer_3d] + else: + self.twod_3D = [self.orbital_viewer_2d] self.on_kpoint_change(None) - - + def _read_arraydata(self, i_spin, i_kpt, i_band): - spin_mult = self.nspin - 1 - kpt_qe_convention = i_kpt + 1 + self.nkpoints_lowres * i_spin*spin_mult - - cube_id = [i for i, f in enumerate(self.list_of_calcs) - if 'K'+str(kpt_qe_convention).zfill(3)+'_'+'B'+str(i_band+1).zfill(3) in list(f)[0]] - - if len(cube_id)==0: + kpt_qe_convention = i_kpt + 1 + self.nkpoints_lowres * i_spin * spin_mult + + cube_id = [ + i + for i, f in enumerate(self.list_of_calcs) + if "K" + + str(kpt_qe_convention).zfill(3) + + "_" + + "B" + + str(i_band + 1).zfill(3) + in list(f)[0] + ] + + if len(cube_id) == 0: return None - - cid=cube_id[0] + + cid = cube_id[0] fname = list(self.list_of_calcs[cid])[0] - - if fname[0] == 'K' and fname[4:6] == '_B': + + if fname[0] == "K" and fname[4:6] == "_B": arraydata = list(self.list_of_calcs[cid])[1] - elif fname.startswith('output_data_multiple'): + elif fname.startswith("output_data_multiple"): arraydata = list(self.list_of_calcs[cid])[1].outputs[fname] else: absfn = list(self.list_of_calcs[cid])[1].outputs.retrieved.open(fname).name with gzip.open(absfn) as fpointer: arraydata = from_cube_to_arraydata(fpointer.read()) - + return arraydata, fname - + def clamp_arraydata_to_zero(self, arraydata): new_ad = copy.deepcopy(arraydata) - data = new_ad.get_array('data') - new_ad.set_array('data', np.maximum(data, 0.0)) + data = new_ad.get_array("data") + new_ad.set_array("data", np.maximum(data, 0.0)) return new_ad - + def on_kpoint_change(self, _=None): """Replot the orbitals in case of the kpoint change.""" - - arraydata_fn = self._read_arraydata(self.bands_viewer.selected_spin, - self.bands_viewer.selected_kpoint-1, - self.bands_viewer.selected_band-1) + + arraydata_fn = self._read_arraydata( + self.bands_viewer.selected_spin, + self.bands_viewer.selected_kpoint - 1, + self.bands_viewer.selected_band - 1, + ) if arraydata_fn is None: with self.output: clear_output() self.info_out.value = "Found no cube files" - #self.orbital_viewer_3d.arraydata = None - #self.orbital_viewer_2d.arraydata = None - display(ipw.VBox([self.info_out], layout=ipw.Layout(margin="200px 0px 0px 0px"))) + # self.orbital_viewer_3d.arraydata = None + # self.orbital_viewer_2d.arraydata = None + display( + ipw.VBox( + [self.info_out], layout=ipw.Layout(margin="200px 0px 0px 0px") + ) + ) else: arraydata, fname = arraydata_fn self.info_out.value = fname - + self.orbital_viewer_2d.export_label = "pk{}_b{}_k{:02}_s{}".format( self._workcalc.pk, self.bands_viewer.selected_band, @@ -786,84 +919,88 @@ def on_kpoint_change(self, _=None): self.bands_viewer.selected_spin, ) self.orbital_viewer_2d.arraydata = self.clamp_arraydata_to_zero(arraydata) - + with self.output: - clear_output() + clear_output() if self.bands_viewer.selected_3D: self.orbital_viewer_3d.arraydata = arraydata hbox = ipw.HBox([self.orbital_viewer_3d]) else: hbox = ipw.HBox([self.orbital_viewer_2d]) - display(ipw.VBox( - [self.info_out, hbox], - layout=ipw.Layout(margin="200px 0px 0px 0px") - )) - + display( + ipw.VBox( + [self.info_out, hbox], + layout=ipw.Layout(margin="200px 0px 0px 0px"), + ) + ) + def gaussian(self, x, fwhm): - sigma = fwhm/2.3548 - return np.exp(-x**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi)) - + sigma = fwhm / 2.3548 + return np.exp(-(x**2) / (2 * sigma**2)) / (sigma * np.sqrt(2 * np.pi)) + def calculate_sts_mapping(self, energy, broadening): - sts_arraydata = None sts_arraydata_meta = None warn = False - - for i_band in range(self.bands_lowres.shape[2]-1, -1, -1): + + for i_band in range(self.bands_lowres.shape[2] - 1, -1, -1): for i_spin in range(self.bands_lowres.shape[0]): for i_kpt in range(self.bands_lowres.shape[1]): bz_w = 1 if i_kpt in [0, 12] else 2 - - orb_energy = self.bands_lowres[i_spin, i_kpt, i_band] - self.vacuum_level - + + orb_energy = ( + self.bands_lowres[i_spin, i_kpt, i_band] - self.vacuum_level + ) + if np.abs(energy - orb_energy) < 2 * broadening: - coef = self.gaussian(energy - orb_energy, broadening) - + ad_out = self._read_arraydata(i_spin, i_kpt, i_band) - + if ad_out is None: warn = True continue - + arraydata, fname = ad_out arraydata = self.clamp_arraydata_to_zero(arraydata) - - orbital_data = arraydata.get_array('data') - + + orbital_data = arraydata.get_array("data") + if sts_arraydata is None: sts_arraydata = bz_w * coef * orbital_data sts_arraydata_meta = copy.deepcopy(arraydata) else: sts_arraydata += bz_w * coef * orbital_data - + if sts_arraydata_meta is None: return None, warn else: - sts_arraydata_meta.set_array('data', sts_arraydata) + sts_arraydata_meta.set_array("data", sts_arraydata) return sts_arraydata_meta, warn - + def on_sts_btn_click(self, _=None): - - sts_arraydata, warn = self.calculate_sts_mapping(self.sts_energy_text.value, self.sts_fwhm_text.value) - + sts_arraydata, warn = self.calculate_sts_mapping( + self.sts_energy_text.value, self.sts_fwhm_text.value + ) + if sts_arraydata is not None: if warn: self.sts_mapping_viewer_wrapper.children = [ ipw.HTML(value="Warning: some relevant bands were not found"), - self.sts_mapping_viewer + self.sts_mapping_viewer, ] else: self.sts_mapping_viewer_wrapper.children = [self.sts_mapping_viewer] - + self.sts_mapping_viewer.export_label = "pk{}_ldos_fwhm{}_e{}".format( self._workcalc.pk, self.sts_fwhm_text.value, self.sts_energy_text.value ) self.sts_mapping_viewer.arraydata = sts_arraydata else: - self.sts_mapping_viewer_wrapper.children = [ipw.HTML(value="Could not find data.")] - - + self.sts_mapping_viewer_wrapper.children = [ + ipw.HTML(value="Could not find data.") + ] + @property def spindensity_calc(self): """Return spindensity plot calculation if present, otherwise return None.""" diff --git a/nanoribbon/viewers/smiles2GNR.py b/nanoribbon/viewers/smiles2GNR.py deleted file mode 100644 index b7961fb..0000000 --- a/nanoribbon/viewers/smiles2GNR.py +++ /dev/null @@ -1,249 +0,0 @@ -# pylint: disable=no-member -"""Widget to convert SMILES to nanoribbons.""" - -import numpy as np -from scipy.stats import mode - -from IPython.display import clear_output -import ipywidgets as ipw -import nglview - -from traitlets import Instance - -from ase import Atoms -from ase.data import covalent_radii,chemical_symbols -from ase.neighborlist import NeighborList -import ase.neighborlist - - -from sklearn import manifold, datasets -from sklearn.decomposition import PCA -import pybel - -class Smiles2GnrWidget(ipw.VBox): - """Conver SMILES into 3D structure.""" - structure = Instance(Atoms, allow_none=True) - SPINNER = """""" - - def __init__(self, title="Smiles to GNR"): - try: - import openbabel # pylint: disable=unused-import - except ImportError: - super().__init__( - [ipw.HTML("The SmilesWidget requires the OpenBabel library, " - "but the library was not found.")]) - return - self.title = title - self.original_structure = None - self.selection = set() - self.smiles = ipw.Text() - - create_structure_btn = ipw.Button(description="Convert SMILES", button_style='info') - create_structure_btn.on_click(self._on_button_pressed) - - self.create_cell_btn = ipw.Button(description="create GNR", button_style='info', disabled=True) - self.create_cell_btn.on_click(self._on_cell_button_pressed) - - self.viewer = nglview.NGLWidget() - self.viewer.stage.set_parameters(mouse_preset='pymol') - self.viewer.observe(self._on_picked, names='picked') - self.select_two = ipw.HTML("") - self.picked_out = ipw.Output() - self.cell_button_out = ipw.Output() - super().__init__([ - self.smiles, - ipw.Label(value="e.g. C1(C2=CC=C(C3=CC=CC=C3)C=C2)=CC=CC=C1"), create_structure_btn, self.select_two, - self.viewer, self.picked_out, self.create_cell_btn, self.cell_button_out - ]) - - @staticmethod - def guess_scaling_factor(atoms): - """Scaling factor to correct the bond length.""" - - # Set bounding box as cell. - c_x = 1.5 * (np.amax(atoms.positions[:, 0]) - np.amin(atoms.positions[:, 0])) - c_y = 1.5 * (np.amax(atoms.positions[:, 1]) - np.amin(atoms.positions[:, 1])) - c_z = 15.0 - atoms.cell = (c_x, c_y, c_z) - atoms.pbc = (True, True, True) - - # Calculate all atom-atom distances. - c_atoms = [a for a in atoms if a.symbol[0] == "C"] - n_atoms = len(c_atoms) - dists = np.zeros([n_atoms, n_atoms]) - for i, atom_a in enumerate(c_atoms): - for j, atom_b in enumerate(c_atoms): - dists[i, j] = np.linalg.norm(atom_a.position - atom_b.position) - - # Find bond distances to closest neighbor. - dists += np.diag([np.inf] * n_atoms) # Don't consider diagonal. - bonds = np.amin(dists, axis=1) - - # Average bond distance. - avg_bond = float(mode(bonds)[0]) - - # Scale box to match equilibrium carbon-carbon bond distance. - cc_eq = 1.4313333333 - return cc_eq / avg_bond - - @staticmethod - def scale(atoms, s): - """Scale atomic positions by the `factor`.""" - c_x, c_y, c_z = atoms.cell - atoms.set_cell((factor * c_x, factor * c_y, c_z), scale_atoms=True) - atoms.center() - return atoms - - @staticmethod - def smiles2d(smiles): - """Create planar molecule from smiles.""" - mol = pybel.readstring("smiles",smiles) - mol.make3D() - #mol.addh() - pybel._builder.Build(mol.OBMol) - #optimize_mol(mol) - f_f = pybel._forcefields["uff"] - if not f_f.Setup(mol.OBMol): - f_f = pybel._forcefields["mmff94"] # pylint: disable=protected-access - if not f_f.Setup(mol.OBMol): - print("Cannot set up forcefield") - return - f_f.Setup(mol.OBMol) - f_f.SteepestDescent(5000, 1.0e-9) - f_f.GetCoordinates(mol.OBMol) - asemol = Atoms() - species=[chemical_symbols[atm.atomicnum] for atm in mol.atoms] - pos=np.asarray([atm.coords for atm in mol.atoms]) - pca = PCA(n_components=3) - pca.fit(pos) - posnew=pca.transform(pos) - posnew[:,2]=0.0 - atoms = Atoms(species, positions=posnew) - return atoms - - @staticmethod - def construct_cell(atoms, id1, id2): - """Construct periodic cell based on two selected equivalent atoms.""" - - pos = [[atoms[id2].x, atoms[id2].y], [atoms[id1].x, atoms[id1].y], [atoms[id2].x, atoms[id1].y]] - - vec = [np.array(pos[0]) - np.array(pos[1]), np.array(pos[2]) - np.array(pos[1])] - c_x = np.linalg.norm(vec[0]) - - angle = np.math.atan2(np.linalg.det([vec[0], vec[1]]), np.dot(vec[0], vec[1])) - if np.abs(angle) > 0.01: - atoms.rotate_euler(center=atoms[id1].position, phi=-angle, theta=0.0, psi=0.0) - - c_y = 15.0 + np.amax(atoms.positions[:, 1]) - np.amin(atoms.positions[:, 1]) - c_z = 15.0 + np.amax(atoms.positions[:, 2]) - np.amin(atoms.positions[:, 2]) - - atoms.cell = (c_x, c_y, c_z) - atoms.pbc = (True, True, True) - atoms.center() - atoms.wrap(eps=0.001) - - # Remove redundant atoms. - tobedel = [] - - n_l = NeighborList([covalent_radii[a.number] for a in atoms], bothways=False, self_interaction=False) - n_l.update(atoms) - - for atm in atoms: - indices, offsets = n_l.get_neighbors(atm.index) - for i, offset in zip(indices, offsets): - dist = np.linalg.norm(atm.position - (atoms.positions[i] + np.dot(offset, atoms.get_cell()))) - if dist < 0.4: - tobedel.append(atoms[i].index) - - del atoms[tobedel] - - # Find unit cell and apply it. - - ## Add Hydrogens. - n_l = NeighborList([covalent_radii[a.number] for a in atoms], bothways=True, self_interaction=False) - n_l.update(atoms) - - need_hydrogen = [] - for atm in atoms: - if len(n_l.get_neighbors(atm.index)[0]) < 3: - if atm.symbol == 'C': - need_hydrogen.append(atm.index) - - print("Added missing Hydrogen atoms: ", need_hydrogen) - - for atm in need_hydrogen: - vec = np.zeros(3) - indices, offsets = n_l.get_neighbors(atoms[atm].index) - for i, offset in zip(indices, offsets): - vec += -atoms[atm].position + (atoms.positions[i] + np.dot(offset, atoms.get_cell())) - vec = -vec / np.linalg.norm(vec) * 1.1 + atoms[atm].position - atoms.append(ase.Atom('H', vec)) - - return atoms - - def _on_picked(self, _=None): - """When an attom is picked.""" - - if 'atom1' not in self.viewer.picked.keys(): - return # did not click on atom - self.create_cell_btn.disabled = True - - with self.picked_out: - clear_output() - - self.viewer.component_0.remove_ball_and_stick() - self.viewer.component_0.remove_ball_and_stick() - self.viewer.add_ball_and_stick() - - idx = self.viewer.picked['atom1']['index'] - - # Toggle. - if idx in self.selection: - self.selection.remove(idx) - else: - self.selection.add(idx) - - if len(self.selection) == 2: - self.create_cell_btn.disabled = False - - #if(selection): - sel_str = ",".join([str(i) for i in sorted(self.selection)]) - print("Selected atoms: " + sel_str) - self.viewer.add_representation('ball+stick', selection="@" + sel_str, color='red', aspectRatio=3.0) - self.viewer.picked = {} # reset, otherwise immidiately selecting same atom again won't create change event - - def _on_button_pressed(self, _=None): - """Convert SMILES to ase structure when button is pressed.""" - - self.create_cell_btn.disabled = True - if not self.smiles.value: - return - - smiles = self.smiles.value.replace(" ", "") - struct = self.smiles2d(smiles) - self.select_two.value = '

Select two equivalent atoms that define the basis vector

' - self.original_structure = struct.copy() - if hasattr(self.viewer, "component_0"): - self.viewer.component_0.remove_ball_and_stick() - cid = self.viewer.component_0.id - self.viewer.remove_component(cid) - - # Empty selection. - self.selection = set() - - # Add new component. - self.viewer.add_component(nglview.ASEStructure(struct)) # adds ball+stick - self.viewer.center() - self.viewer.handle_resize() - - def _on_cell_button_pressed(self, _=None): - """Generate GNR button pressed.""" - with self.cell_button_out: - clear_output() - if len(self.selection) != 2: - print("You must select exactly two atoms") - return - - id1 = sorted(self.selection)[0] - id2 = sorted(self.selection)[1] - self.structure = self.construct_cell(self.original_structure, id1, id2) \ No newline at end of file diff --git a/nanoribbon/viewers/smiles2gnr.py b/nanoribbon/viewers/smiles2gnr.py deleted file mode 100644 index b7961fb..0000000 --- a/nanoribbon/viewers/smiles2gnr.py +++ /dev/null @@ -1,249 +0,0 @@ -# pylint: disable=no-member -"""Widget to convert SMILES to nanoribbons.""" - -import numpy as np -from scipy.stats import mode - -from IPython.display import clear_output -import ipywidgets as ipw -import nglview - -from traitlets import Instance - -from ase import Atoms -from ase.data import covalent_radii,chemical_symbols -from ase.neighborlist import NeighborList -import ase.neighborlist - - -from sklearn import manifold, datasets -from sklearn.decomposition import PCA -import pybel - -class Smiles2GnrWidget(ipw.VBox): - """Conver SMILES into 3D structure.""" - structure = Instance(Atoms, allow_none=True) - SPINNER = """""" - - def __init__(self, title="Smiles to GNR"): - try: - import openbabel # pylint: disable=unused-import - except ImportError: - super().__init__( - [ipw.HTML("The SmilesWidget requires the OpenBabel library, " - "but the library was not found.")]) - return - self.title = title - self.original_structure = None - self.selection = set() - self.smiles = ipw.Text() - - create_structure_btn = ipw.Button(description="Convert SMILES", button_style='info') - create_structure_btn.on_click(self._on_button_pressed) - - self.create_cell_btn = ipw.Button(description="create GNR", button_style='info', disabled=True) - self.create_cell_btn.on_click(self._on_cell_button_pressed) - - self.viewer = nglview.NGLWidget() - self.viewer.stage.set_parameters(mouse_preset='pymol') - self.viewer.observe(self._on_picked, names='picked') - self.select_two = ipw.HTML("") - self.picked_out = ipw.Output() - self.cell_button_out = ipw.Output() - super().__init__([ - self.smiles, - ipw.Label(value="e.g. C1(C2=CC=C(C3=CC=CC=C3)C=C2)=CC=CC=C1"), create_structure_btn, self.select_two, - self.viewer, self.picked_out, self.create_cell_btn, self.cell_button_out - ]) - - @staticmethod - def guess_scaling_factor(atoms): - """Scaling factor to correct the bond length.""" - - # Set bounding box as cell. - c_x = 1.5 * (np.amax(atoms.positions[:, 0]) - np.amin(atoms.positions[:, 0])) - c_y = 1.5 * (np.amax(atoms.positions[:, 1]) - np.amin(atoms.positions[:, 1])) - c_z = 15.0 - atoms.cell = (c_x, c_y, c_z) - atoms.pbc = (True, True, True) - - # Calculate all atom-atom distances. - c_atoms = [a for a in atoms if a.symbol[0] == "C"] - n_atoms = len(c_atoms) - dists = np.zeros([n_atoms, n_atoms]) - for i, atom_a in enumerate(c_atoms): - for j, atom_b in enumerate(c_atoms): - dists[i, j] = np.linalg.norm(atom_a.position - atom_b.position) - - # Find bond distances to closest neighbor. - dists += np.diag([np.inf] * n_atoms) # Don't consider diagonal. - bonds = np.amin(dists, axis=1) - - # Average bond distance. - avg_bond = float(mode(bonds)[0]) - - # Scale box to match equilibrium carbon-carbon bond distance. - cc_eq = 1.4313333333 - return cc_eq / avg_bond - - @staticmethod - def scale(atoms, s): - """Scale atomic positions by the `factor`.""" - c_x, c_y, c_z = atoms.cell - atoms.set_cell((factor * c_x, factor * c_y, c_z), scale_atoms=True) - atoms.center() - return atoms - - @staticmethod - def smiles2d(smiles): - """Create planar molecule from smiles.""" - mol = pybel.readstring("smiles",smiles) - mol.make3D() - #mol.addh() - pybel._builder.Build(mol.OBMol) - #optimize_mol(mol) - f_f = pybel._forcefields["uff"] - if not f_f.Setup(mol.OBMol): - f_f = pybel._forcefields["mmff94"] # pylint: disable=protected-access - if not f_f.Setup(mol.OBMol): - print("Cannot set up forcefield") - return - f_f.Setup(mol.OBMol) - f_f.SteepestDescent(5000, 1.0e-9) - f_f.GetCoordinates(mol.OBMol) - asemol = Atoms() - species=[chemical_symbols[atm.atomicnum] for atm in mol.atoms] - pos=np.asarray([atm.coords for atm in mol.atoms]) - pca = PCA(n_components=3) - pca.fit(pos) - posnew=pca.transform(pos) - posnew[:,2]=0.0 - atoms = Atoms(species, positions=posnew) - return atoms - - @staticmethod - def construct_cell(atoms, id1, id2): - """Construct periodic cell based on two selected equivalent atoms.""" - - pos = [[atoms[id2].x, atoms[id2].y], [atoms[id1].x, atoms[id1].y], [atoms[id2].x, atoms[id1].y]] - - vec = [np.array(pos[0]) - np.array(pos[1]), np.array(pos[2]) - np.array(pos[1])] - c_x = np.linalg.norm(vec[0]) - - angle = np.math.atan2(np.linalg.det([vec[0], vec[1]]), np.dot(vec[0], vec[1])) - if np.abs(angle) > 0.01: - atoms.rotate_euler(center=atoms[id1].position, phi=-angle, theta=0.0, psi=0.0) - - c_y = 15.0 + np.amax(atoms.positions[:, 1]) - np.amin(atoms.positions[:, 1]) - c_z = 15.0 + np.amax(atoms.positions[:, 2]) - np.amin(atoms.positions[:, 2]) - - atoms.cell = (c_x, c_y, c_z) - atoms.pbc = (True, True, True) - atoms.center() - atoms.wrap(eps=0.001) - - # Remove redundant atoms. - tobedel = [] - - n_l = NeighborList([covalent_radii[a.number] for a in atoms], bothways=False, self_interaction=False) - n_l.update(atoms) - - for atm in atoms: - indices, offsets = n_l.get_neighbors(atm.index) - for i, offset in zip(indices, offsets): - dist = np.linalg.norm(atm.position - (atoms.positions[i] + np.dot(offset, atoms.get_cell()))) - if dist < 0.4: - tobedel.append(atoms[i].index) - - del atoms[tobedel] - - # Find unit cell and apply it. - - ## Add Hydrogens. - n_l = NeighborList([covalent_radii[a.number] for a in atoms], bothways=True, self_interaction=False) - n_l.update(atoms) - - need_hydrogen = [] - for atm in atoms: - if len(n_l.get_neighbors(atm.index)[0]) < 3: - if atm.symbol == 'C': - need_hydrogen.append(atm.index) - - print("Added missing Hydrogen atoms: ", need_hydrogen) - - for atm in need_hydrogen: - vec = np.zeros(3) - indices, offsets = n_l.get_neighbors(atoms[atm].index) - for i, offset in zip(indices, offsets): - vec += -atoms[atm].position + (atoms.positions[i] + np.dot(offset, atoms.get_cell())) - vec = -vec / np.linalg.norm(vec) * 1.1 + atoms[atm].position - atoms.append(ase.Atom('H', vec)) - - return atoms - - def _on_picked(self, _=None): - """When an attom is picked.""" - - if 'atom1' not in self.viewer.picked.keys(): - return # did not click on atom - self.create_cell_btn.disabled = True - - with self.picked_out: - clear_output() - - self.viewer.component_0.remove_ball_and_stick() - self.viewer.component_0.remove_ball_and_stick() - self.viewer.add_ball_and_stick() - - idx = self.viewer.picked['atom1']['index'] - - # Toggle. - if idx in self.selection: - self.selection.remove(idx) - else: - self.selection.add(idx) - - if len(self.selection) == 2: - self.create_cell_btn.disabled = False - - #if(selection): - sel_str = ",".join([str(i) for i in sorted(self.selection)]) - print("Selected atoms: " + sel_str) - self.viewer.add_representation('ball+stick', selection="@" + sel_str, color='red', aspectRatio=3.0) - self.viewer.picked = {} # reset, otherwise immidiately selecting same atom again won't create change event - - def _on_button_pressed(self, _=None): - """Convert SMILES to ase structure when button is pressed.""" - - self.create_cell_btn.disabled = True - if not self.smiles.value: - return - - smiles = self.smiles.value.replace(" ", "") - struct = self.smiles2d(smiles) - self.select_two.value = '

Select two equivalent atoms that define the basis vector

' - self.original_structure = struct.copy() - if hasattr(self.viewer, "component_0"): - self.viewer.component_0.remove_ball_and_stick() - cid = self.viewer.component_0.id - self.viewer.remove_component(cid) - - # Empty selection. - self.selection = set() - - # Add new component. - self.viewer.add_component(nglview.ASEStructure(struct)) # adds ball+stick - self.viewer.center() - self.viewer.handle_resize() - - def _on_cell_button_pressed(self, _=None): - """Generate GNR button pressed.""" - with self.cell_button_out: - clear_output() - if len(self.selection) != 2: - print("You must select exactly two atoms") - return - - id1 = sorted(self.selection)[0] - id2 = sorted(self.selection)[1] - self.structure = self.construct_cell(self.original_structure, id1, id2) \ No newline at end of file diff --git a/nanoribbon/viewers/structures.py b/nanoribbon/viewers/structures.py deleted file mode 100644 index c2ec8e7..0000000 --- a/nanoribbon/viewers/structures.py +++ /dev/null @@ -1,685 +0,0 @@ -"""Module to provide functionality to import structures.""" - -import os -import tempfile -import datetime -from collections import OrderedDict -from traitlets import Bool -import ipywidgets as ipw - -from aiida.orm import CalcFunctionNode, CalcJobNode, Node, QueryBuilder, WorkChainNode, StructureData -from .utils import get_ase_from_file - - -class StructureManagerWidget(ipw.VBox): # pylint: disable=too-many-instance-attributes - '''Upload a structure and store it in AiiDA database. - - Useful class members: - :ivar has_structure: whether the widget contains a structure - :vartype has_structure: bool - :ivar frozen: whenter the widget is frozen (can't be modified) or not - :vartype frozen: bool - :ivar structure_node: link to AiiDA structure object - :vartype structure_node: StructureData or CifData''' - - has_structure = Bool(False) - frozen = Bool(False) - DATA_FORMATS = ('StructureData', 'CifData') - - def __init__(self, importers, storable=True, node_class=None, **kwargs): - """ - :param storable: Whether to provide Store button (together with Store format) - :type storable: bool - :param node_class: AiiDA node class for storing the structure. - Possible values: 'StructureData', 'CifData' or None (let the user decide). - Note: If your workflows require a specific node class, better fix it here. - :param examples: list of tuples each containing a name and a path to an example structure - :type examples: list - :param importers: list of tuples each containing a name and an object for data importing. Each object - should containt an empty `on_structure_selection()` method that has two parameters: structure_ase, name - :type examples: list""" - - from .viewers import StructureDataViewer - if not importers: # we make sure the list is not empty - raise ValueError("The parameter importers should contain a list (or tuple) of tuples " - "(\"importer name\", importer), got a falsy object.") - - self.structure_ase = None - self._structure_node = None - - self.viewer = StructureDataViewer(downloadable=False) - - self.btn_store = ipw.Button(description='Store in AiiDA', disabled=True) - self.btn_store.on_click(self._on_click_store) - - # Description that will is stored along with the new structure. - self.structure_description = ipw.Text(placeholder="Description (optional)") - - # Select format to store in the AiiDA database. - self.data_format = ipw.RadioButtons(options=self.DATA_FORMATS, description='Data type:') - self.data_format.observe(self.reset_structure, names=['value']) - - if len(importers) == 1: - # If there is only one importer - no need to make tabs. - self._structure_sources_tab = importers[0][1] - # Assigning a function which will be called when importer provides a structure. - importers[0][1].on_structure_selection = self.select_structure - else: - self._structure_sources_tab = ipw.Tab() # Tabs. - self._structure_sources_tab.children = [i[1] for i in importers] # One importer per tab. - for i, (label, importer) in enumerate(importers): - # Labeling tabs. - self._structure_sources_tab.set_title(i, label) - # Assigning a function which will be called when importer provides a structure. - importer.on_structure_selection = self.select_structure - - if storable: - if node_class is None: - store = [self.btn_store, self.data_format, self.structure_description] - elif node_class not in self.DATA_FORMATS: - raise ValueError("Unknown data format '{}'. Options: {}".format(node_class, self.DATA_FORMATS)) - else: - self.data_format.value = node_class - store = [self.btn_store, self.structure_description] - else: - store = [self.structure_description] - store = ipw.HBox(store) - - super().__init__(children=[self._structure_sources_tab, self.viewer, store], **kwargs) - - def reset_structure(self, change=None): # pylint: disable=unused-argument - if self.frozen: - return - self._structure_node = None - self.viewer.structure = None - - def select_structure(self, structure_ase, name): - """Select structure - - :param structure_ase: ASE object containing structure - :type structure_ase: ASE Atoms - :param name: File name with extension but without path - :type name: str""" - - if self.frozen: - return - self._structure_node = None - if not structure_ase: - self.btn_store.disabled = True - self.has_structure = False - self.structure_ase = None - self.structure_description.value = '' - self.reset_structure() - return - self.btn_store.disabled = False - self.has_structure = True - self.structure_description.value = "{} ({})".format(structure_ase.get_chemical_formula(), name) - self.structure_ase = structure_ase - self.viewer.structure = structure_ase - - def _on_click_store(self, change): # pylint: disable=unused-argument - self.store_structure() - - def store_structure(self, label=None, description=None): - """Stores the structure in AiiDA database.""" - if self.frozen: - return - if self.structure_node is None: - return - if self.structure_node.is_stored: - print("Already stored in AiiDA: " + repr(self.structure_node) + " skipping..") - return - if label: - self.structure_node.label = label - if description: - self.structure_node.description = description - self.structure_node.store() - print("Stored in AiiDA: " + repr(self.structure_node)) - - def freeze(self): - """Do not allow any further modifications""" - self._structure_sources_tab.layout.visibility = 'hidden' - self.frozen = True - self.btn_store.disabled = True - self.structure_description.disabled = True - self.data_format.disabled = True - - @property - def node_class(self): - return self.data_format.value - - @node_class.setter - def node_class(self, value): - if self.frozen: - return - self.data_format.value = value - - @property - def structure_node(self): - """Returns AiiDA StructureData node.""" - if self._structure_node is None: - if self.structure_ase is None: - return None - # perform conversion - if self.data_format.value == 'CifData': - from aiida.orm.nodes.data.cif import CifData - self._structure_node = CifData() - self._structure_node.set_ase(self.structure_ase) - else: # Target format is StructureData - self._structure_node = StructureData(ase=self.structure_ase) - self._structure_node.description = self.structure_description.value - self._structure_node.label = self.structure_ase.get_chemical_formula() - return self._structure_node - - -class StructureUploadWidget(ipw.VBox): - """Class that allows to upload structures from user's computer.""" - - def __init__(self, text="Upload Structure"): - from fileupload import FileUploadWidget - - self.on_structure_selection = lambda structure_ase, name: None - self.file_path = None - self.file_upload = FileUploadWidget(text) - supported_formats = ipw.HTML( - """ - Supported structure formats - """) - self.file_upload.observe(self._on_file_upload, names='data') - super().__init__(children=[self.file_upload, supported_formats]) - - def _on_file_upload(self, change): # pylint: disable=unused-argument - """When file upload button is pressed.""" - self.file_path = os.path.join(tempfile.mkdtemp(), self.file_upload.filename) - with open(self.file_path, 'w') as fobj: - fobj.write(self.file_upload.data.decode("utf-8")) - structure_ase = get_ase_from_file(self.file_path) - self.on_structure_selection(structure_ase=structure_ase, name=self.file_upload.filename) - - -class StructureExamplesWidget(ipw.VBox): - """Class to provide example structures for selection.""" - - def __init__(self, examples, **kwargs): - self.on_structure_selection = lambda structure_ase, name: None - self._select_structure = ipw.Dropdown(options=self.get_example_structures(examples)) - self._select_structure.observe(self._on_select_structure, names=['value']) - super().__init__(children=[self._select_structure], **kwargs) - - @staticmethod - def get_example_structures(examples): - """Get the list of example structures.""" - if not isinstance(examples, list): - raise ValueError("parameter examples should be of type list, {} given".format(type(examples))) - return [("Select structure", False)] + examples - - def _on_select_structure(self, change): # pylint: disable=unused-argument - """When structure is selected.""" - if not self._select_structure.value: - return - structure_ase = get_ase_from_file(self._select_structure.value) - self.on_structure_selection(structure_ase=structure_ase, name=self._select_structure.label) - - -class StructureBrowserWidget(ipw.VBox): - """Class to query for structures stored in the AiiDA database.""" - - def __init__(self): - # Find all process labels - qbuilder = QueryBuilder() - qbuilder.append(WorkChainNode, project="label") - qbuilder.order_by({WorkChainNode: {'ctime': 'desc'}}) - process_labels = {i[0] for i in qbuilder.all() if i[0]} - - layout = ipw.Layout(width="900px") - self.mode = ipw.RadioButtons(options=['all', 'uploaded', 'edited', 'calculated'], - layout=ipw.Layout(width="25%")) - - # Date range - self.dt_now = datetime.datetime.now() - self.dt_end = self.dt_now - datetime.timedelta(days=10) - self.date_start = ipw.Text(value='', description='From: ', style={'description_width': '120px'}) - - self.date_end = ipw.Text(value='', description='To: ') - self.date_text = ipw.HTML(value='

Select the date range:

') - self.btn_date = ipw.Button(description='Search', layout={'margin': '1em 0 0 0'}) - self.age_selection = ipw.VBox( - [self.date_text, ipw.HBox([self.date_start, self.date_end]), self.btn_date], - layout={ - 'border': '1px solid #fafafa', - 'padding': '1em' - }) - - # Labels - self.drop_label = ipw.Dropdown(options=({'All'}.union(process_labels)), - value='All', - description='Process Label', - style={'description_width': '120px'}, - layout={'width': '50%'}) - - self.btn_date.on_click(self.search) - self.mode.observe(self.search, names='value') - self.drop_label.observe(self.search, names='value') - - h_line = ipw.HTML('
') - box = ipw.VBox([self.age_selection, h_line, ipw.HBox([self.mode, self.drop_label])]) - - self.results = ipw.Dropdown(layout=layout) - self.results.observe(self._on_select_structure) - self.search() - super(StructureBrowserWidget, self).__init__([box, h_line, self.results]) - - @staticmethod - def preprocess(): - """Search structures in AiiDA database.""" - queryb = QueryBuilder() - queryb.append(StructureData, filters={'extras': {'!has_key': 'formula'}}) - for itm in queryb.all(): # iterall() would interfere with set_extra() - formula = itm[0].get_formula() - itm[0].set_extra("formula", formula) - - def search(self, change=None): # pylint: disable=unused-argument - """Launch the search of structures in AiiDA database.""" - self.preprocess() - - qbuild = QueryBuilder() - try: # If the date range is valid, use it for the search - self.start_date = datetime.datetime.strptime(self.date_start.value, '%Y-%m-%d') - self.end_date = datetime.datetime.strptime(self.date_end.value, '%Y-%m-%d') + datetime.timedelta(hours=24) - except ValueError: # Otherwise revert to the standard (i.e. last 7 days) - self.start_date = self.dt_end - self.end_date = self.dt_now + datetime.timedelta(hours=24) - - self.date_start.value = self.start_date.strftime('%Y-%m-%d') - self.date_end.value = self.end_date.strftime('%Y-%m-%d') - - filters = {} - filters['ctime'] = {'and': [{'<=': self.end_date}, {'>': self.start_date}]} - if self.drop_label.value != 'All': - qbuild.append(WorkChainNode, filters={'label': self.drop_label.value}) - # print(qbuild.all()) - # qbuild.append(CalcJobNode, with_incoming=WorkChainNode) - qbuild.append(StructureData, with_incoming=WorkChainNode, filters=filters) - else: - if self.mode.value == "uploaded": - qbuild2 = QueryBuilder() - qbuild2.append(StructureData, project=["id"]) - qbuild2.append(Node, with_outgoing=StructureData) - processed_nodes = [n[0] for n in qbuild2.all()] - if processed_nodes: - filters['id'] = {"!in": processed_nodes} - qbuild.append(StructureData, filters=filters) - - elif self.mode.value == "calculated": - qbuild.append(CalcJobNode) - qbuild.append(StructureData, with_incoming=CalcJobNode, filters=filters) - - elif self.mode.value == "edited": - qbuild.append(CalcFunctionNode) - qbuild.append(StructureData, with_incoming=CalcFunctionNode, filters=filters) - - elif self.mode.value == "all": - qbuild.append(StructureData, filters=filters) - - qbuild.order_by({StructureData: {'ctime': 'desc'}}) - matches = {n[0] for n in qbuild.iterall()} - matches = sorted(matches, reverse=True, key=lambda n: n.ctime) - - options = OrderedDict() - options["Select a Structure ({} found)".format(len(matches))] = False - - for mch in matches: - label = "PK: %d" % mch.pk - label += " | " + mch.ctime.strftime("%Y-%m-%d %H:%M") - label += " | " + mch.get_extra("formula") - label += " | " + mch.description - options[label] = mch - - self.results.options = options - - def _on_select_structure(self, change): # pylint: disable=unused-argument - """When a structure was selected.""" - if not self.results.value: - return - structure_ase = self.results.value.get_ase() - formula = structure_ase.get_chemical_formula() - if self.on_structure_selection is not None: - self.on_structure_selection(structure_ase=structure_ase, name=formula) - - def on_structure_selection(self, structure_ase, name): - pass - - -class SmilesWidget(ipw.VBox): - """Conver SMILES into 3D structure.""" - - SPINNER = """""" - - def __init__(self): - try: - import openbabel # pylint: disable=unused-import - except ImportError: - super().__init__( - [ipw.HTML("The SmilesWidget requires the OpenBabel library, " - "but the library was not found.")]) - return - - self.smiles = ipw.Text() - self.create_structure_btn = ipw.Button(description="Generate molecule", button_style='info') - self.create_structure_btn.on_click(self._on_button_pressed) - self.output = ipw.HTML("") - super().__init__([self.smiles, self.create_structure_btn, self.output]) - - @staticmethod - def pymol_2_ase(pymol): - """Convert pymol object into ASE Atoms.""" - import numpy as np - from ase import Atoms, Atom - from ase.data import chemical_symbols - - asemol = Atoms() - for atm in pymol.atoms: - asemol.append(Atom(chemical_symbols[atm.atomicnum], atm.coords)) - asemol.cell = np.amax(asemol.positions, axis=0) - np.amin(asemol.positions, axis=0) + [10] * 3 - asemol.pbc = True - asemol.center() - return asemol - - def _optimize_mol(self, mol): - """Optimize a molecule using force field (needed for complex SMILES).""" - - # Note, the pybel module imported below comes together with openbabel package. Do not confuse it with - # pybel package available on PyPi: https://pypi.org/project/pybel/ - import pybel # pylint:disable=import-error - - self.output.value = "Screening possible conformers {}".format(self.SPINNER) #font-size:20em; - - f_f = pybel._forcefields["mmff94"] # pylint: disable=protected-access - if not f_f.Setup(mol.OBMol): - f_f = pybel._forcefields["uff"] # pylint: disable=protected-access - if not f_f.Setup(mol.OBMol): - self.output.value = "Cannot set up forcefield" - return - - # initial cleanup before the weighted search - f_f.SteepestDescent(5500, 1.0e-9) - f_f.WeightedRotorSearch(15000, 500) - f_f.ConjugateGradients(6500, 1.0e-10) - f_f.GetCoordinates(mol.OBMol) - self.output.value = "" - - def _on_button_pressed(self, change): # pylint: disable=unused-argument - """Convert SMILES to ase structure when button is pressed.""" - self.output.value = "" - - # Note, the pybel module imported below comes together with openbabel package. Do not confuse it with - # pybel package available on PyPi: https://pypi.org/project/pybel/ - import pybel # pylint:disable=import-error - - if not self.smiles.value: - return - - mol = pybel.readstring("smi", self.smiles.value) - self.output.value = """SMILES to 3D conversion {}""".format(self.SPINNER) - mol.make3D() - - pybel._builder.Build(mol.OBMol) # pylint: disable=protected-access - mol.addh() - self._optimize_mol(mol) - - structure_ase = self.pymol_2_ase(mol) - formula = structure_ase.get_chemical_formula() - if self.on_structure_selection is not None: - self.on_structure_selection(structure_ase=structure_ase, name=formula) - - def on_structure_selection(self, structure_ase, name): - pass - -import numpy as np -from scipy.stats import mode -from numpy.linalg import norm -from pysmiles import read_smiles,write_smiles -from rdkit.Chem.rdmolfiles import MolFromSmiles,MolToMolFile -import networkx as nx -import math -from ase import Atoms -from ase.visualize import view -from IPython.display import display, clear_output -import ipywidgets as ipw -import nglview -from ase.data import covalent_radii -from ase.neighborlist import NeighborList -import ase.neighborlist - -class SmilesWidget(ipw.VBox): - """Conver SMILES into 3D structure.""" - - SPINNER = """""" - - def __init__(self): - try: - import openbabel # pylint: disable=unused-import - except ImportError: - super().__init__( - [ipw.HTML("The SmilesWidget requires the OpenBabel library, " - "but the library was not found.")]) - return - - self.selection = set() - self.cell_ready = False - self.smiles = ipw.Text() - self.create_structure_btn = ipw.Button(description="Convert SMILES", button_style='info') - self.create_structure_btn.on_click(self._on_button_pressed) - self.create_cell_btn = ipw.Button(description="create GNR", button_style='info') - self.create_cell_btn.on_click(self._on_button2_pressed) - self.viewer = nglview.NGLWidget() - self.viewer.observe(self._on_picked, names='picked') - self.output = ipw.HTML("") - self.picked_out = ipw.Output() - self.button2_out = ipw.Output() - super().__init__([self.smiles, self.create_structure_btn,self.viewer,self_picked_out, self.output,self.button2_out]) - -######## - @staticmethod - def guess_scaling_factor(atoms): - import numpy as np - # set bounding box as cell - cx = 1.5 * (np.amax(atoms.positions[:,0]) - np.amin(atoms.positions[:,0])) - cy = 1.5 * (np.amax(atoms.positions[:,1]) - np.amin(atoms.positions[:,1])) - cz = 15.0 - atoms.cell = (cx, cy, cz) - atoms.pbc = (True,True,True) - - # calculate all atom-atom distances - c_atoms = [a for a in atoms if a.symbol[0]=="C"] - n = len(c_atoms) - dists = np.zeros([n,n]) - for i, a in enumerate(c_atoms): - for j, b in enumerate(c_atoms): - dists[i,j] = norm(a.position - b.position) - - # find bond distances to closest neighbor - dists += np.diag([np.inf]*n) # don't consider diagonal - bonds = np.amin(dists, axis=1) - - # average bond distance - avg_bond = float(mode(bonds)[0]) - - # scale box to match equilibrium carbon-carbon bond distance - cc_eq = 1.4313333333 - s = cc_eq / avg_bond - return s - - @staticmethod - def scale(atoms, s): - cx, cy, cz = atoms.cell - atoms.set_cell((s*cx, s*cy, cz), scale_atoms=True) - atoms.center() - return atoms - - @staticmethod - def smiles2D(smiles): - mol = MolFromSmiles(smiles) - from rdkit.Chem import AllChem - # generate the 2D coordinates - AllChem.Compute2DCoords(mol) - # get the 2D coordinates - for c in mol.GetConformers(): - coords=c.GetPositions() - # get the atom labels - ll=[] - for i in mol.GetAtoms(): - #ll.append(i.GetSymbol()) - ll.append(i.GetAtomicNum()) - ll=np.asarray(ll) - # create an ASE frame - c=Atoms('{:d}N'.format(len(coords))) - c.set_positions(coords) - c.set_atomic_numbers(ll) - return c - - @staticmethod - def construct_cell(atoms, id1, id2): - - p1 = [atoms[id1].x, atoms[id1].y] - p0 = [atoms[id2].x, atoms[id2].y] - p2 = [atoms[id2].x, atoms[id1].y] - - v0 = np.array(p0) - np.array(p1) - v1 = np.array(p2) - np.array(p1) - - angle = np.math.atan2(np.linalg.det([v0,v1]),np.dot(v0,v1)) - - #angle=np.degrees(angle) - - cx = norm(v0) - - #print np.degrees(angle),v0,v1,p0,p1 - if np.abs(angle) > 0.01: - # s.euler_rotate(phi=angle,theta=0,psi=0,center(x[id1],y[id1],z[id1])) - atoms.rotate_euler(center=atoms[id1].position, phi=-angle,theta=0.0,psi=0.0) - - yrange = np.amax(atoms.positions[:,1])-np.amin(atoms.positions[:,1]) - zrange = np.amax(atoms.positions[:,2])-np.amin(atoms.positions[:,2]) - cy = 15.0 + yrange - cz = 15.0 + zrange - - atoms.cell = (cx,cy,cz) - atoms.pbc = (True,True,True) - atoms.center() - atoms.wrap(eps=0.001) - - #### REMOVE REDUNDANT ATOMS - tobedel = [] - - cov_radii = [covalent_radii[a.number] for a in atoms] - nl = NeighborList(cov_radii, bothways = False, self_interaction = False) - nl.update(atoms) - - for a in atoms: - indices, offsets = nl.get_neighbors(a.index) - for i, offset in zip(indices, offsets): - dist = norm(a.position -(atoms.positions[i] + np.dot(offset, atoms.get_cell()))) - if dist < 0.4 : - tobedel.append(atoms[i].index) - - del atoms[tobedel] - - #### ENDFIND UNIT CELL AND APPLIES IT - - #### ADD Hydrogens - cov_radii = [covalent_radii[a.number] for a in atoms] - nl = NeighborList(cov_radii, bothways = True, self_interaction = False) - nl.update(atoms) - - need_a_H = [] - for a in atoms: - nlist=nl.get_neighbors(a.index)[0] - if len(nlist)<3: - if a.symbol=='C': - need_a_H.append(a.index) - - print("Added missing Hydrogen atoms: ", need_a_H) - - dCH=1.1 - for a in need_a_H: - vec = np.zeros(3) - indices, offsets = nl.get_neighbors(atoms[a].index) - for i, offset in zip(indices, offsets): - vec += -atoms[a].position +(atoms.positions[i] + np.dot(offset, atoms.get_cell())) - vec = -vec/norm(vec)*dCH - vec += atoms[a].position - htoadd = ase.Atom('H',vec) - atoms.append(htoadd) - - return atoms - - def _on_picked(self,ca): - - self.cell_ready = False - - if 'atom1' not in self.viewer.picked.keys(): - return # did not click on atom - with self.picked_out: - clear_output() - - #viewer.clear_representations() - self.viewer.component_0.remove_ball_and_stick() - self.viewer.component_0.remove_ball_and_stick() - self.viewer.add_ball_and_stick() - #viewer.add_unitcell() - - idx = self.viewer.picked['atom1']['index'] - - # toggle - if idx in self.selection: - self.selection.remove(idx) - else: - self.selection.add(idx) - - #if(selection): - sel_str = ",".join([str(i) for i in sorted(self.selection)]) - print("Selected atoms: "+ sel_str) - self.viewer.add_representation('ball+stick', selection="@"+sel_str, color='red', aspectRatio=3.0) - #else: - # print ("nothing selected") - self.viewer.picked = {} # reset, otherwise immidiately selecting same atom again won't create change event - - - - def _on_button_pressed(self, change): # pylint: disable=unused-argument - """Convert SMILES to ase structure when button is pressed.""" - self.output.value = "" - - - if not self.smiles.value: - return - - smiles=self.smiles.value.replace(" ", "") - c=self.smiles2D(smiles) - # set the cell - - scaling_fac=self.guess_scaling_factor(c) - scaled_structure=self.scale(c,scaling_fac) - self.original_structure=c.copy() - - def _on_button2_pressed(self, change): - with self.button2_out: - clear_output() - self.cell_ready = False - - if len(self.selection) != 2: - print("You must select exactly two atoms") - return - - - id1 = sorted(self.selection)[0] - id2 = sorted(self.selection)[1] - new_structure = self.construct_cell(self.original_structure, id1, id2) - formula = new_structure.get_chemical_formula() - if self.on_structure_selection is not None: - self.on_structure_selection(structure_ase=new_structure, name=formula) - self.cell_ready = True - - def on_structure_selection(self, structure_ase, name): - pass diff --git a/nanoribbon/viewers/utils.py b/nanoribbon/viewers/utils.py index 54489bb..cb66791 100644 --- a/nanoribbon/viewers/utils.py +++ b/nanoribbon/viewers/utils.py @@ -1,17 +1,17 @@ """Utility functions for the nanoribbon workchain viewers.""" -import numpy as np -import matplotlib.pyplot as plt -import matplotlib.colors as mc import colorsys -from ase.data import covalent_radii, atomic_numbers -from ase.neighborlist import NeighborList -from ase.data.colors import cpk_colors +import matplotlib.colors as mc +import matplotlib.pyplot as plt +import numpy as np +from aiida.orm import CalcJobNode, QueryBuilder, WorkChainNode # AiiDA imports from aiida.plugins import DataFactory -from aiida.orm import CalcJobNode, QueryBuilder, WorkChainNode +from ase.data import atomic_numbers, covalent_radii +from ase.data.colors import cpk_colors +from ase.neighborlist import NeighborList # AiiDA data types. ArrayData = DataFactory("array") # pylint: disable=invalid-name @@ -26,8 +26,8 @@ def get_calc_by_label(workcalc, label): def get_calcs_by_label(workcalc, label): """Get step calculation of a workchain by its name.""" qbld = QueryBuilder() - qbld.append(WorkChainNode, filters={'uuid': workcalc.uuid}) - qbld.append(CalcJobNode, with_incoming=WorkChainNode, filters={'label': label}) + qbld.append(WorkChainNode, filters={"uuid": workcalc.uuid}) + qbld.append(CalcJobNode, with_incoming=WorkChainNode, filters={"label": label}) calcs = [c[0] for c in qbld.all()] for calc in calcs: assert calc.is_finished_ok @@ -38,7 +38,9 @@ def from_cube_to_arraydata(cube_content): """Convert cube file to the AiiDA ArrayData object.""" lines = cube_content.splitlines() natoms = int(lines[2].split()[0]) # The number of atoms listed in the file - header = lines[:6 + natoms] # Header of the file: comments, the voxel, and the number of atoms and datapoints + header = lines[ + : 6 + natoms + ] # Header of the file: comments, the voxel, and the number of atoms and datapoints # Parse the declared dimensions of the volumetric data x_line = header[3].split() @@ -50,8 +52,13 @@ def from_cube_to_arraydata(cube_content): # Get the vectors describing the basis voxel voxel_array = np.array( - [[x_line[1], x_line[2], x_line[3]], [y_line[1], y_line[2], y_line[3]], [z_line[1], z_line[2], z_line[3]]], - dtype=np.float64) + [ + [x_line[1], x_line[2], x_line[3]], + [y_line[1], y_line[2], y_line[3]], + [z_line[1], z_line[2], z_line[3]], + ], + dtype=np.float64, + ) atm_numbers = np.empty(natoms, int) coordinates = np.empty((natoms, 3)) for i in range(natoms): @@ -62,21 +69,22 @@ def from_cube_to_arraydata(cube_content): # Get the volumetric data data_array = np.empty(xdim * ydim * zdim, dtype=float) cursor = 0 - for line in lines[6 + natoms:]: # The actual data: atoms and volumetric data + for line in lines[6 + natoms :]: # The actual data: atoms and volumetric data lsplitted = line.split() - data_array[cursor:cursor + len(lsplitted)] = lsplitted + data_array[cursor : cursor + len(lsplitted)] = lsplitted cursor += len(lsplitted) arraydata = ArrayData() - arraydata.set_array('voxel', voxel_array) - arraydata.set_array('data', data_array.reshape((xdim, ydim, zdim))) - arraydata.set_array('data_units', np.array('e/bohr^3')) - arraydata.set_array('coordinates_units', np.array('bohr')) - arraydata.set_array('coordinates', coordinates) - arraydata.set_array('atomic_numbers', atm_numbers) + arraydata.set_array("voxel", voxel_array) + arraydata.set_array("data", data_array.reshape((xdim, ydim, zdim))) + arraydata.set_array("data_units", np.array("e/bohr^3")) + arraydata.set_array("coordinates_units", np.array("bohr")) + arraydata.set_array("coordinates", coordinates) + arraydata.set_array("atomic_numbers", atm_numbers) return arraydata + def adjust_lightness(color, amount=0.9): try: c = mc.cnames[color] @@ -93,7 +101,7 @@ def plot_struct_2d(ax_plt, atoms, alpha): # Plot overlayed structure. strct = atoms.repeat((4, 1, 1)) - strct.positions[:, 0] -= atoms.cell[0,0] + strct.positions[:, 0] -= atoms.cell[0, 0] cov_radii = [covalent_radii[a.number] for a in strct] nlist = NeighborList(cov_radii, bothways=True, self_interaction=False) nlist.update(strct) @@ -103,21 +111,27 @@ def plot_struct_2d(ax_plt, atoms, alpha): pos = atm.position nmbrs = atomic_numbers[atm.symbol] ax_plt.add_artist( - plt.Circle((pos[0], pos[1]), - covalent_radii[nmbrs] * 0.4, - color=adjust_lightness(cpk_colors[nmbrs], amount=0.90), - fill=True, - clip_on=True, - alpha=alpha)) + plt.Circle( + (pos[0], pos[1]), + covalent_radii[nmbrs] * 0.4, + color=adjust_lightness(cpk_colors[nmbrs], amount=0.90), + fill=True, + clip_on=True, + alpha=alpha, + ) + ) # Bonds. for theneig in nlist.get_neighbors(atm.index)[0]: pos = (strct[theneig].position + atm.position) / 2 pos0 = atm.position - if (pos[0] - pos0[0])**2 + (pos[1] - pos0[1])**2 < 2: - ax_plt.plot([pos0[0], pos[0]], [pos0[1], pos[1]], - color=adjust_lightness(cpk_colors[nmbrs], amount=0.90), - linewidth=2, - linestyle='-', - solid_capstyle='butt', - alpha=alpha) + if (pos[0] - pos0[0]) ** 2 + (pos[1] - pos0[1]) ** 2 < 2: + ax_plt.plot( + [pos0[0], pos[0]], + [pos0[1], pos[1]], + color=adjust_lightness(cpk_colors[nmbrs], amount=0.90), + linewidth=2, + linestyle="-", + solid_capstyle="butt", + alpha=alpha, + ) diff --git a/nanoribbon/workchains/__init__.py b/nanoribbon/workchains/__init__.py deleted file mode 100644 index 2b2f800..0000000 --- a/nanoribbon/workchains/__init__.py +++ /dev/null @@ -1,7 +0,0 @@ -"""Nanoribbon work chains.""" -# pylint: disable=unused-import - -from __future__ import absolute_import - -from .nanoribbonwork import NanoribbonWorkChain -from .ks_orbitals import KSWorkChain diff --git a/nanoribbon/workchains/ks_orbitals.py b/nanoribbon/workchains/ks_orbitals.py deleted file mode 100644 index 8aec130..0000000 --- a/nanoribbon/workchains/ks_orbitals.py +++ /dev/null @@ -1,255 +0,0 @@ -import six -import numpy as np -import itertools - -# AiiDA imports -from aiida.orm import Code, Computer, Dict, Int, Float, KpointsData, Str, StructureData, SinglefileData, load_node -from aiida.engine import WorkChain, ToContext, CalcJob, run, submit -#from aiida.orm.nodes.data.upf import get_pseudos_dict, get_pseudos_from_structure - -# aiida_quantumespresso imports -from aiida.engine import ExitCode -from aiida_quantumespresso.calculations.pw import PwCalculation -from aiida_quantumespresso.calculations.pp import PpCalculation -from aiida_quantumespresso.calculations.projwfc import ProjwfcCalculation -from aiida_quantumespresso.utils.pseudopotential import validate_and_prepare_pseudos_inputs - -# aditional imports -from . import aux_script_strings - - -class KSWorkChain(WorkChain): - - @classmethod - def define(cls, spec): - super(KSWorkChain, cls).define(spec) - spec.input("pw_code", valid_type=Code) - spec.input("pp_code", valid_type=Code) - spec.input("workchain", valid_type=Int) - spec.input("inpkpoints", valid_type=Str) - spec.input("inpbands", valid_type=Str) - # TODO: check why it does not work - #spec.inputs("metadata.label", valid_type=six.string_types, - # default="NanoribbonWorkChain", non_db=True, help="Label of the work chain.") - spec.outline( - cls.run_scf, - cls.run_bands, - cls.run_export_orbitals - ) - #spec.dynamic_output() - spec.outputs.dynamic = True - - - def run_scf(self): - self.ctx.orig_w = load_node(self.inputs.workchain.value) - self.ctx.pseudo_family=self.ctx.orig_w.inputs.pseudo_family - self.ctx.structure = self.ctx.orig_w.called_descendants[-2].outputs.output_structure - kpoints = self.ctx.orig_w.called_descendants[-3].inputs.kpoints - nkpt=kpoints.get_kpoints_mesh()[0][0] - return self._submit_pw_calc(self.ctx.structure, label="scf", runtype='scf', - kpoints=kpoints,nkpt=nkpt, wallhours=1) - - - def run_bands(self): - prev_calc = self.ctx.scf - self._check_prev_calc(prev_calc) - parent_folder = prev_calc.outputs.remote_folder - kpoints=self.ctx.orig_w.called_descendants[-7].inputs.kpoints - return self._submit_pw_calc(self.ctx.structure, - label="bands", - parent_folder=parent_folder, - runtype='bands', - kpoints=kpoints, - nkpt=12, - wallhours=1) - - - # ========================================================================= - def run_export_orbitals(self): - self.report("Running pp.x to export KS orbitals") - builder = PpCalculation.get_builder() - builder.code = self.inputs.pp_code - nproc_mach = builder.code.computer.get_default_mpiprocs_per_machine() - - prev_calc = self.ctx.bands - self._check_prev_calc(prev_calc) - builder.parent_folder = prev_calc.outputs.remote_folder - - nel = prev_calc.res.number_of_electrons - bands=list(map(int, self.inputs.inpbands.value.split())) - kpoints=list(map(int, self.inputs.inpkpoints.value.split())) - #nkpt = prev_calc.res.number_of_k_points - #nbnd = prev_calc.res.number_of_bands - #nspin = prev_calc.res.number_of_spin_components - volume = prev_calc.res.volume - nhours = int(1) - - nnodes=int(prev_calc.attributes['resources']['num_machines']) - npools = int(prev_calc.inputs.settings['cmdline'][1]) - #nproc_mach=int(prev_calc.attributes['resources']['num_mpiprocs_per_machine']) - for ib,ik in itertools.product(bands,kpoints): - builder.parameters = Dict(dict={ - 'inputpp': { - # contribution of a selected wavefunction - # to charge density - 'plot_num': 7, - 'kpoint(1)': ik, - 'kpoint(2)': ik, - 'kband(1)': ib+int(nel/2), - 'kband(2)': ib+int(nel/2), - }, - 'plot': { - 'iflag': 3, # 3D plot - 'output_format': 6, # CUBE format - 'fileout': '_orbital.cube', - }, - }) - - builder.metadata.label = "export_orbitals" - builder.metadata.options = { - "resources": { - "num_machines": nnodes, - "num_mpiprocs_per_machine": nproc_mach, - }, - "max_wallclock_seconds": nhours * 60 * 60, # 6 hours - "append_text": aux_script_strings.cube_cutter, - "withmpi": True, - } - - - builder.settings = Dict( - dict={'additional_retrieve_list': ['*.cube.gz'], - 'cmdline': - ["-npools", str(npools)] - } - ) - running = self.submit(builder) - label = 'export_orbitals_{}_{}'.format(ib,ik) - self.to_context(**{label:running}) - return - - - # ========================================================================= - def _check_prev_calc(self, prev_calc): - error = None - output_fname = prev_calc.attributes['output_filename'] - if not prev_calc.is_finished_ok: - error = "Previous calculation failed" #in state: "+prev_calc.get_state() - elif output_fname not in prev_calc.outputs.retrieved.list_object_names(): - error = "Previous calculation did not retrive {}".format(output_fname) - else: - content = prev_calc.outputs.retrieved.get_object_content(output_fname) - if "JOB DONE." not in content: - error = "Previous calculation not DONE." - if error: - self.report("ERROR: "+error) - #self.abort(msg=error) ## ABORT WILL NOT WORK, not defined - #raise Exception(error) - return ExitCode(450) - - # ========================================================================= - def _submit_pw_calc(self, structure, label, runtype, - kpoints=None,nkpt=None, wallhours=24, parent_folder=None): - self.report("Running pw.x for "+label) - builder = PwCalculation.get_builder() - - builder.code = self.inputs.pw_code - builder.structure = structure - builder.parameters = self._get_parameters(structure, runtype,label) - builder.pseudos = validate_and_prepare_pseudos_inputs(structure, None, self.ctx.pseudo_family) - - - if parent_folder: - builder.parent_folder = parent_folder - - # kpoints - use_symmetry = runtype != "bands" - #kpoints = self._get_kpoints(nkpoints, use_symmetry=use_symmetry) - builder.kpoints = kpoints - - # parallelization settings - ## TEMPORARY double pools in case of spin - spinpools=int(1) - start_mag = self._get_magnetization(structure) - if any([m != 0 for m in start_mag.values()]): - spinpools = int(2) - npools = spinpools*min( int(nkpt/5), int(5) ) - natoms = len(structure.sites) - nnodes = (1 + int(natoms/60) ) * npools - - builder.metadata.label = label - builder.metadata.options = { - "resources": { "num_machines": nnodes }, - "withmpi": True, - "max_wallclock_seconds": wallhours * 60 * 60, - } - - builder.settings = Dict(dict={'cmdline': ["-npools", str(npools)]}) - - future = self.submit(builder) - return ToContext(**{label:future}) - - # ========================================================================= - def _get_parameters(self, structure, runtype, label): - params = {'CONTROL': { - 'calculation': runtype, - 'wf_collect': True, - 'forc_conv_thr': 0.0001, - 'nstep': 500, - }, - 'SYSTEM': { - 'ecutwfc': 50., - 'ecutrho': 400., - 'occupations': 'smearing', - 'degauss': 0.001, - }, - 'ELECTRONS': { - 'conv_thr': 1.e-8, - 'mixing_beta': 0.25, - 'electron_maxstep': 50, - 'scf_must_converge': False, - }, - } - - if label == 'cell_opt1': - params['CONTROL']['forc_conv_thr']=0.0005 - if runtype == "vc-relax": - # in y and z direction there is only vacuum - params['CELL'] = {'cell_dofree': 'x'} - - # if runtype == "bands": - # params['CONTROL']['restart_mode'] = 'restart' - - start_mag = self._get_magnetization(structure) - if any([m != 0 for m in start_mag.values()]): - params['SYSTEM']['nspin'] = 2 - params['SYSTEM']['starting_magnetization'] = start_mag - - return Dict(dict=params) - - # ========================================================================= - def _get_kpoints(self, nx, use_symmetry=True): - nx = max(1, nx) - - kpoints = KpointsData() - if use_symmetry: - kpoints.set_kpoints_mesh([nx, 1, 1], offset=[0.0, 0.0, 0.0]) - else: - # list kpoints explicitly - points = [[r, 0.0, 0.0] for r in np.linspace(0, 0.5, nx)] - kpoints.set_kpoints(points) - - return kpoints - - - # ========================================================================= - def _get_magnetization(self, structure): - start_mag = {} - for i in structure.kinds: - if i.name.endswith("1"): - start_mag[i.name] = 1.0 - elif i.name.endswith("2"): - start_mag[i.name] = -1.0 - else: - start_mag[i.name] = 0.0 - return start_mag diff --git a/search.ipynb b/search.ipynb index 967cd15..6de33da 100644 --- a/search.ipynb +++ b/search.ipynb @@ -19,9 +19,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "display(NanoribbonSearchWidget())" diff --git a/setup.cfg b/setup.cfg index 0a49ca1..f29f4d0 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,38 +1,47 @@ -[aiidalab] -title = Empa nanotech@surfaces Laboratory - Graphene nanoribbons - [metadata] name = nanoribbons version = 1.5.1 -author = nanotech@surfaces laboratory, Empa description = Tools for graphene nanoribbons, developed at the nanotech@surfaces laboratory, Empa long_description = file: README.md long_description_content_type = text/markdown url = https://github.com/nanotech-empa/aiidalab-empa-nanoribbons -project_urls = - Logo = https://raw.githubusercontent.com/nanotech-empa/aiidalab-empa-nanoribbons/develop/miscellaneous/logos/empa.png - Bug Tracker = https://github.com/nanotech-empa/aiidalab-empa-nanoribbons/issues +author = nanotech@surfaces laboratory, Empa +license = MIT +license_file = LICENSE.txt classifiers = - Programming Language :: Python :: 3 + Development Status :: 5 - Production/Stable License :: OSI Approved :: MIT License Operating System :: OS Independent - Development Status :: 5 - Production/Stable + Programming Language :: Python :: 3 + Programming Language :: Python :: 3 :: Only +project_urls = + Logo = https://raw.githubusercontent.com/nanotech-empa/aiidalab-empa-nanoribbons/develop/miscellaneous/logos/empa.png + Bug Tracker = https://github.com/nanotech-empa/aiidalab-empa-nanoribbons/issues [options] packages = find: -python_requires = >=3.7 install_requires = aiida-core~=1.0 - aiida-quantumespresso~=3.5 aiida-nanotech-empa~=0.6 + aiida-quantumespresso~=3.5 aiidalab-widgets-base~=1.0 ase +python_requires = >=3.7 [options.extras_require] dev = bumpver==2021.1118 pre-commit==2.11.1 +[aiidalab] +title = Empa nanotech@surfaces Laboratory - Graphene nanoribbons + +[flake8] +ignore = + E501 + W503 + E203 + [bumpver] current_version = "v1.5.1" version_pattern = "vMAJOR.MINOR.PATCH[PYTAGNUM]" diff --git a/setup.py b/setup.py index c1057cf..b908cbe 100644 --- a/setup.py +++ b/setup.py @@ -1,3 +1,3 @@ import setuptools -setuptools.setup() \ No newline at end of file +setuptools.setup() diff --git a/start.py b/start.py index 6deb5e0..c7e03b3 100644 --- a/start.py +++ b/start.py @@ -1,21 +1,23 @@ import ipywidgets as ipw + def get_start_widget(appbase, jupbase): - #http://fontawesome.io/icons/ + # http://fontawesome.io/icons/ template = """ - + - +
""" - + html = template.format(appbase=appbase, jupbase=jupbase) return ipw.HTML(html) - -#EOF + + +# EOF diff --git a/submit_ks.ipynb b/submit_ks.ipynb deleted file mode 100644 index f7c0cd4..0000000 --- a/submit_ks.ipynb +++ /dev/null @@ -1,124 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Submit additional KS for nanoribbon" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%aiida" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# General imports\n", - "import numpy as np\n", - "import ipywidgets as ipw\n", - "from collections import OrderedDict\n", - "from ase.data import vdw_radii\n", - "from IPython.display import display, clear_output, HTML\n", - "import itertools\n", - "\n", - "# AiiDA & AiiDA lab imports\n", - "from aiida.orm import SinglefileData\n", - "from aiidalab_widgets_base import CodeDropdown, StructureManagerWidget, StructureBrowserWidget, StructureUploadWidget, SubmitButtonWidget, SmilesWidget\n", - "from aiida.engine import submit\n", - "\n", - "# Work Chains\n", - "KSWorkChain = WorkflowFactory('ksnanoribbon')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "style = {'description_width': '120px'}\n", - "layout = {'width': '70%'}\n", - "original_workchain = ipw.Text(description='PK',style={'description_width': '70px'}, layout={'width': '14%'})\n", - "kpoints = ipw.Text(description='k-points: 1-12 U; 13-24 D',style={'description_width': '70px'}, layout={'width': '14%'})\n", - "bands = ipw.Text(description='bands: 0 for nel/2',style={'description_width': '70px'}, layout={'width': '14%'})\n", - "\n", - "pw_code_dropdown = CodeDropdown(input_plugin='quantumespresso.pw')\n", - "pp_code_dropdown = CodeDropdown(input_plugin='quantumespresso.pp')\n", - "\n", - "\n", - "display(ipw.VBox([original_workchain,kpoints,bands,pw_code_dropdown,pp_code_dropdown]))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def on_submit():\n", - " with submit_out:\n", - " clear_output()\n", - " if not pw_code_dropdown.selected_code:\n", - " print(\"Please select a pw code.\")\n", - " return\n", - " if not pp_code_dropdown.selected_code:\n", - " print(\"Please select a pp code.\")\n", - " return\n", - "\n", - " #pw_code, pp_code, projwfc_code = qe_code_groups[drop_codes.value]\n", - " builder = KSWorkChain.get_builder()\n", - " builder.pw_code = pw_code_dropdown.selected_code\n", - " builder.pp_code = pp_code_dropdown.selected_code\n", - " builder.inpkpoints = Str(kpoints.value)\n", - " builder.inpbands = Str(bands.value)\n", - " builder.workchain = Int(original_workchain.value)\n", - " builder.metadata = {\n", - " \"description\": 'additional KS',\n", - " \"label\": \"KSWorkChain\",\n", - " }\n", - " return builder" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "btn_submit = SubmitButtonWidget(process=KSWorkChain,widgets_values=on_submit)\n", - "submit_out = ipw.Output()\n", - "display(btn_submit, submit_out)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.4" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} From 71c25e480394c1f6e1dd10e37f86d5eda4d3aebd Mon Sep 17 00:00:00 2001 From: Aliaksandr Yakutovich Date: Wed, 31 May 2023 09:12:57 +0000 Subject: [PATCH 2/4] Fix almost all pre-commit errors. --- nanoribbon/editors/__init__.py | 4 - nanoribbon/viewers/__init__.py | 5 - nanoribbon/viewers/cdxml2gnr.py | 73 +-------------- nanoribbon/viewers/igor.py | 64 +++++++------ nanoribbon/viewers/pdos_computed.py | 34 +++---- nanoribbon/viewers/search.py | 136 +++++++++++++++++----------- nanoribbon/viewers/show_computed.py | 31 +++---- nanoribbon/viewers/utils.py | 2 +- 8 files changed, 148 insertions(+), 201 deletions(-) diff --git a/nanoribbon/editors/__init__.py b/nanoribbon/editors/__init__.py index 5cd1e9f..c1aac6c 100644 --- a/nanoribbon/editors/__init__.py +++ b/nanoribbon/editors/__init__.py @@ -1,9 +1,5 @@ """Nanoribbon AiiDAlab viewers.""" -from aiida import load_profile - -load_profile() - from .replicate import NanoribbonReplicateEditor __all__ = ["NanoribbonReplicateEditor"] diff --git a/nanoribbon/viewers/__init__.py b/nanoribbon/viewers/__init__.py index 4b36fbf..8109910 100644 --- a/nanoribbon/viewers/__init__.py +++ b/nanoribbon/viewers/__init__.py @@ -1,10 +1,5 @@ """Nanoribbon AiiDAlab viewers.""" - -from aiida import load_profile - -load_profile() - from .cdxml2gnr import CdxmlUpload2GnrWidget from .pdos_computed import NanoribbonPDOSWidget from .search import NanoribbonSearchWidget diff --git a/nanoribbon/viewers/cdxml2gnr.py b/nanoribbon/viewers/cdxml2gnr.py index 97349af..586fee7 100644 --- a/nanoribbon/viewers/cdxml2gnr.py +++ b/nanoribbon/viewers/cdxml2gnr.py @@ -28,7 +28,7 @@ class CdxmlUpload2GnrWidget(ipw.VBox): def __init__(self, title="CDXML to GNR", description="Upload Structure"): try: - import openbabel # pylint: disable=unused-import + import openbabel # noqa: F401 except ImportError: super().__init__( [ @@ -165,73 +165,6 @@ def construct_cell(atoms, id1, id2): atoms.center() atoms.wrap(eps=0.001) - # #### REMOVE EXTRA ATOMS - # some_too_close=True - # while some_too_close and len(atoms)>1: - # #ite+=1 - # some_too_close=False - # cov_radii = [covalent_radii[a.number] for a in atoms] - # nl = NeighborList(cov_radii, bothways = True, self_interaction = False) - # nl.update(atoms) - # for a in atoms: - # if some_too_close: - # break - # - # indices, offsets = nl.get_neighbors(a.index) - # #print(ite,a,indices) - # for i, offset in zip(indices, offsets): - # dist = np.linalg.norm(a.position -(atoms.positions[i] + np.dot(offset, atoms.get_cell()))) - # if dist < 0.4 : - # some_too_close=True - # #remove H atoms before others - # if atoms[i].symbol != 'H' and a.symbol == 'H': - # del atoms[a.index] - # #print('deleted ',a.index) - # else: - # del atoms[i] - # #print('deleted ',i) - # break - # - # elif a.symbol == "C" : #this is executed only if dist > 0.4 - # - # #this part removes H from overcoordinated C - # #if atoms[i].symbol=='H' and [ida.symbol for ida in atoms[indices]].count('C') == 3: - # - # #this part removes H from all C and N - # if atoms[i].symbol=='H' or a.symbol=='N' : - # some_too_close=True - # del atoms[i] - # #print('deleted ',i) - # break - # - # - # #### ADD Hydrogens - # cov_radii = [covalent_radii[a.number] for a in atoms] - # nl = NeighborList(cov_radii, bothways = True, self_interaction = False) - # nl.update(atoms) - # - # need_a_H = [] - # for a in atoms: - # nlist=nl.get_neighbors(a.index)[0] - # if len(nlist)<3: - # if a.symbol=='C' or a.symbol=='N': - # need_a_H.append(a.index) - # - # #print("Added missing Hydrogen atoms: ", need_a_H) - # - # dCH=1.1 - # for a in need_a_H: - # vec = np.zeros(3) - # indices, offsets = nl.get_neighbors(atoms[a].index) - # for i, offset in zip(indices, offsets): - # vec += -atoms[a].position +(atoms.positions[i] + np.dot(offset, atoms.get_cell())) - # vec = -vec/np.linalg.norm(vec)*dCH - # vec += atoms[a].position - # htoadd = ase.Atom('H',vec) - # atoms.append(htoadd) - # - # return atoms - # Remove redundant atoms. ORIGINAL tobedel = [] @@ -256,7 +189,7 @@ def construct_cell(atoms, id1, id2): # Find unit cell and apply it. - ## Add Hydrogens. + # Add Hydrogens. n_l = NeighborList( [covalent_radii[a.number] for a in atoms], bothways=True, @@ -339,7 +272,7 @@ def _on_file_upload(self, change=None): self.mols[molid] = m listmols.append( (str(molid) + ": " + m.formula, molid) - ) ## m MUST BE a pb object!!! + ) # m must be a pb object!!! molid += 1 self.allmols.options = listmols self.allmols.value = 0 diff --git a/nanoribbon/viewers/igor.py b/nanoribbon/viewers/igor.py index d50f623..066a552 100644 --- a/nanoribbon/viewers/igor.py +++ b/nanoribbon/viewers/igor.py @@ -9,12 +9,22 @@ import numpy as np +class FileNotBeginWithIgorError(OSError): + def __init__(self): + super().__init__("File does not begin with 'IGOR'") + + +class MissingBeginStatementError(OSError): + def __init__(self): + super().__init__("Missing 'BEGIN' statement of data block") + + class Axis: """Represents an axis of an IGOR wave""" - def __init__(self, symbol, min, delta, unit, wavename=None): + def __init__(self, symbol, minimum, delta, unit, wavename=None): self.symbol = symbol - self.min = min + self.min = minimum self.delta = delta self.unit = unit self.wavename = wavename @@ -90,7 +100,7 @@ def read(self, fname): line = lines.pop(0) if not line == "IGOR": - raise OSError("Files does not begin with 'IGOR'") + raise FileNotBeginWithIgorError() line = lines.pop(0) while not re.match("WAVES", line): @@ -102,9 +112,9 @@ def read(self, fname): line = lines.pop(0) if not line == "BEGIN": - raise OSError("Missing 'BEGIN' statement of data block") + raise MissingBeginStatementError() - # read data + # Read data. datastring = "" line = lines.pop(0) while not re.match("END", line): @@ -113,7 +123,7 @@ def read(self, fname): data = np.array(datastring.split(), dtype=float) self.data = data.reshape(grid) - # read axes + # Read axes. line = lines.pop(0) matches = re.findall("SetScale.+?(?:;|$)", line) self.axes = [] @@ -122,10 +132,6 @@ def read(self, fname): ax.read(match) self.axes.append(ax) - # the rest is discarded... - # line = lines.pop(0) - # print(line) - @property def extent(self): """Returns extent for plotting""" @@ -152,12 +158,12 @@ def write(self, fname): class Wave1d(Wave): """1d Igor wave""" - default_parameters = dict( - xmin=0.0, - xdelta=None, - xlabel="x", - ylabel="y", - ) + default_parameters = { + "xmin": 0.0, + "xdelta": None, + "xlabel": "x", + "ylabel": "y", + } def __init__(self, data=None, axes=None, name="1d", **kwargs): """Initialize 1d IGOR wave""" @@ -174,7 +180,7 @@ def __init__(self, data=None, axes=None, name="1d", **kwargs): p = self.parameters x = Axis( symbol="x", - min=p["xmin"], + minimum=p["xmin"], delta=p["xdelta"], unit=p["xlabel"], wavename=self.name, @@ -192,16 +198,16 @@ def print_data(self): class Wave2d(Wave): """2d Igor wave""" - default_parameters = dict( - xmin=0.0, - xdelta=None, - xmax=None, - xlabel="x", - ymin=0.0, - ydelta=None, - ymax=None, - ylabel="y", - ) + default_parameters = { + "xmin": 0.0, + "xdelta": None, + "xmax": None, + "xlabel": "x", + "ymin": 0.0, + "ydelta": None, + "ymax": None, + "ylabel": "y", + } def __init__(self, data=None, axes=None, name=None, **kwargs): """Initialize 2d Igor wave @@ -238,14 +244,14 @@ def __init__(self, data=None, axes=None, name=None, **kwargs): x = Axis( symbol="x", - min=p["xmin"], + minimum=p["xmin"], delta=p["xdelta"], unit=p["xlabel"], wavename=self.name, ) y = Axis( symbol="y", - min=p["ymin"], + minimum=p["ymin"], delta=p["ydelta"], unit=p["ylabel"], wavename=self.name, diff --git a/nanoribbon/viewers/pdos_computed.py b/nanoribbon/viewers/pdos_computed.py index 4d8fb79..8b35e5e 100644 --- a/nanoribbon/viewers/pdos_computed.py +++ b/nanoribbon/viewers/pdos_computed.py @@ -11,7 +11,7 @@ import scipy.ndimage from aiida.orm import CalcJobNode, QueryBuilder, WorkChainNode from IPython.core.display import HTML -from IPython.display import clear_output +from IPython.display import clear_output, display from matplotlib.ticker import FormatStrFormatter on_band_click_global = None @@ -29,7 +29,7 @@ def get_calcs_by_label(workcalc, label): qb.append(CalcJobNode, with_incoming=WorkChainNode, filters={"label": label}) calcs = [c[0] for c in qb.all()] for calc in calcs: - assert calc.is_finished_ok == True + assert calc.is_finished_ok return calcs @@ -128,8 +128,6 @@ def __init__(self, workcalc, **kwargs): self.bands = np.swapaxes(self.eigvalues, 1, 2) + self.vacuum_level - ####BOH - style = {"description_width": "200px"} layout = ipw.Layout(width="600px") self.sigma_slider = ipw.FloatSlider( @@ -285,10 +283,10 @@ def parse_old_atomic_proj_xml(self, root): ) for i in range(self.nspins): for k in range(self.nkpoints): - for l in range(self.natwfcs): + for j in range(self.natwfcs): spintag = "SPIN.%d/" % (i + 1) if self.nspins > 1 else "" raw = root.find( - "PROJECTIONS/K-POINT.%d/%sATMWFC.%d" % (k + 1, spintag, l + 1) + "PROJECTIONS/K-POINT.%d/%sATMWFC.%d" % (k + 1, spintag, j + 1) ).text arr = np.fromstring(raw.replace(",", "\n"), sep="\n") arr2 = arr.reshape( @@ -297,13 +295,13 @@ def parse_old_atomic_proj_xml(self, root): arr3 = np.sum( np.square(arr2), axis=1 ) # calculate square of abs value - self.projections[i, :, k, l] = arr3 + self.projections[i, :, k, j] = arr3 # 4 - def calc_pdos(self, sigma, ngauss, Emin, Emax, atmwfcs=None): - DeltaE = 0.01 - x = np.arange(Emin, Emax, DeltaE) + def calc_pdos(self, sigma, ngauss, emin, emax, atmwfcs=None): + delta_e = 0.01 + x = np.arange(emin, emax, delta_e) # calculate histogram for all spins, bands, and kpoints in parallel xx = np.tile( @@ -326,7 +324,7 @@ def calc_pdos(self, sigma, ngauss, Emin, Emax, atmwfcs=None): # 5 def igor_pdos(self): center = (self.homo + self.lumo) / 2.0 - Emin, Emax = center - 3.0, center + 3.0 + emin, emax = center - 3.0, center + 3.0 if self.selected_atoms: atmwfcs = [ k - 1 @@ -338,8 +336,8 @@ def igor_pdos(self): pdos = self.calc_pdos( ngauss=self.ngauss_slider.value, sigma=self.sigma_slider.value, - Emin=Emin, - Emax=Emax, + emin=emin, + emax=emax, atmwfcs=atmwfcs, ) e = pdos[0] @@ -438,14 +436,12 @@ def plot_pdos(self, ax, pdos_full, ispin, pdos=None): ax.set_xlabel("DOS [a.u.]") ax.xaxis.set_major_formatter(FormatStrFormatter("%.0f")) - if pdos != None: + if pdos is not None: x, y = pdos col = matplotlib.colors.to_rgb(self.colorpicker.value) ax.plot(y[:, ispin], x, color="k") - # ax.plot(y[:,ispin], x, color='blue') tfrm = matplotlib.transforms.Affine2D().rotate_deg(90) + ax.transData ax.fill_between(x, 0.0, -y[:, ispin], facecolor=col, transform=tfrm) - # ax.fill_between(x, 0.0, -y[:,ispin], facecolor='cyan', transform=tfrm) # 7 @@ -465,7 +461,7 @@ def plot_bands(self, ax, ispin, fig_aspect, atmwfcs=None): y_data = y_datas[:, i_band] ax.plot(x_data, y_data, "-", color="black") - ### plot the projection on bands + # Plot the projection on bands. if atmwfcs is not None: line_widths = np.zeros(len(x_data)) for atomwfc in atmwfcs: @@ -514,7 +510,7 @@ def plot_all(self): fig_aspect = figsize[1] / (figsize[0] / 4.0) * 0.5 / (emax - emin) sharey = None - pdos_full = self.calc_pdos(ngauss=ngauss, sigma=sigma, Emin=emin, Emax=emax) + pdos_full = self.calc_pdos(ngauss=ngauss, sigma=sigma, emin=emin, emax=emax) # DOS projected to selected atoms pdos = None @@ -528,7 +524,7 @@ def plot_all(self): ] print("Selected atmwfcs: " + str(atmwfcs)) pdos = self.calc_pdos( - ngauss=ngauss, sigma=sigma, Emin=emin, Emax=emax, atmwfcs=atmwfcs + ngauss=ngauss, sigma=sigma, emin=emin, emax=emax, atmwfcs=atmwfcs ) for ispin in range(self.nspins): diff --git a/nanoribbon/viewers/search.py b/nanoribbon/viewers/search.py index 9faca0f..ecc5e2b 100644 --- a/nanoribbon/viewers/search.py +++ b/nanoribbon/viewers/search.py @@ -13,6 +13,60 @@ HA2EV = 27.211386245988 +class WrongCodeError(Exception): + """Raised when facing an unexpected code behavior.""" + + def __init__(self): + super().__init__("Something wrong has been implemented. Revise the code!") + + +class FermiEnergyAndBandsEnergiesError(ValueError): + """Raised when the Fermi energy is below or above all the bands energies.""" + + def __init__(self, where): + super().__init__( + "The Fermi energy is {where} the bands energies. I don't know what to do." + ) + + +class NeedMoreBandsError(ValueError): + def __init__(self): + super().__init__( + "To understand if it is a metal or insulator, I need more bands than number of electrons." + ) + + +class FermiEnergyOrOccupationsNotPresentError(KeyError): + def __init__(self): + super().__init__( + "I cannot determine metallicity, since I don't have either fermi energy or occupations." + ) + + +class EitherNumberOfElectronsOrFermiEnergyError(ValueError): + def __init__(self): + super().__init__( + "Specify either the number of electrons or the Fermi energy, but not both." + ) + + +class FilesNotFoundError(Exception): + def __init__(self): + super().__init__( + "Did not find 'vacuum_hartree.dat' in the repository or 'output_data' array in the output." + ) + + +class NoAiidaOutError(Exception): + def __init__(self, label): + super().__init__(f"Calculation {label} did not retrive aiida.out") + + +class CalculationNotFinishedCorrectlyError(Exception): + def __init__(self, label, state): + super().__init__(f"Calculation {label} in state {state}.") + + class NanoribbonSearchWidget(ipw.VBox): STYLE = {"description_width": "120px"} LAYOUT = ipw.Layout(width="80%") @@ -38,12 +92,12 @@ def __init__(self, **kwargs): style=self.STYLE, ) - def slider(desc, min, max): + def slider(desc, mininum, maximum): return ipw.FloatRangeSlider( description=desc, - min=min, - max=max, - value=[min, max], + min=mininum, + max=maximum, + value=[mininum, maximum], step=0.01, layout=self.LAYOUT, style=self.STYLE, @@ -174,23 +228,20 @@ def get_calc_by_label(workcalc, label): for label in all_steps: calc = get_calc_by_label(workcalc, label) if calc.process_state.value != "finished": - raise ( - Exception( - f"Calculation {label} in state {calc.process_state.value}." - ) + raise CalculationNotFinishedCorrectlyError( + label=label, state=calc.process_state.value ) if calc.attributes["output_filename"] not in [ obj.name for obj in calc.outputs.retrieved.list_objects() ]: - raise (Exception(f"Calculation {label} did not retrive aiida.out")) + raise NoAiidaOutError(label=label) content = calc.outputs.retrieved.get_object_content( calc.attributes["output_filename"] ) if "JOB DONE." not in content: - # raise(Exception("Calculation {} did not print JOB DONE.".format(label))) print(f"Calculation {label} did not print JOB DONE.") # energies @@ -219,7 +270,6 @@ def get_calc_by_label(workcalc, label): # HOMO, LUMO, and Gap bands_calc = get_calc_by_label(workcalc, "bands") bands = bands_calc.outputs.output_band - parts = self.find_bandgap(bands, fermi_energy=fermi_energy) is_insulator, gap, homo, lumo = self.find_bandgap( bands, fermi_energy=fermi_energy ) @@ -238,15 +288,13 @@ def get_calc_by_label(workcalc, label): except FileNotFoundError: try: data = export_hartree_calc.outputs.output_data.get_array("data") - except KeyError: - raise Exception( - "Did not find 'vacuum_hartree.dat' file in the file repository or" - "'output_data' array in the output." - ) + except KeyError as exc: + raise FilesNotFoundError() from exc + vacuum_level = np.mean(data) * HA2EV * 0.5 workcalc.set_extra("vacuum_level", vacuum_level) - # store shifted energies + # Store shifted energies. workcalc.set_extra("fermi_energy", fermi_energy - vacuum_level) if is_insulator: workcalc.set_extra("homo", homo - vacuum_level) @@ -345,10 +393,7 @@ def nint(num): return int(num - 0.5) if fermi_energy and number_electrons: - raise ValueError( - "Specify either the number of electrons or the " - "Fermi energy, but not both" - ) + raise EitherNumberOfElectronsOrFermiEnergyError() assert bandsdata.units == "eV" stored_bands = bandsdata.get_bands() @@ -358,7 +403,7 @@ def nint(num): # spin up and spin down array # put all spins on one band per kpoint - bands = np.concatenate([_ for _ in stored_bands], axis=1) + bands = np.concatenate(list(stored_bands), axis=1) else: bands = stored_bands @@ -369,11 +414,8 @@ def nint(num): if number_electrons is None: try: _, stored_occupations = bandsdata.get_bands(also_occupations=True) - except KeyError: - raise KeyError( - "Cannot determine metallicity if I don't have " - "either fermi energy, or occupations" - ) + except KeyError as exc: + raise FermiEnergyOrOccupationsNotPresentError() from exc # put the occupations in the same order of bands, also in case of multiple bands if len(stored_occupations.shape) == 3: @@ -381,9 +423,7 @@ def nint(num): # spin up and spin down array # put all spins on one band per kpoint - occupations = np.concatenate( - [_ for _ in stored_occupations], axis=1 - ) + occupations = np.concatenate(list(stored_occupations), axis=1) else: occupations = stored_occupations @@ -423,11 +463,8 @@ def nint(num): homo = [_[0][_[1]] for _ in zip(bands, homo_indexes)] try: lumo = [_[0][_[1] + 1] for _ in zip(bands, homo_indexes)] - except IndexError: - raise ValueError( - "To understand if it is a metal or insulator, " - "need more bands than n_band=number_electrons" - ) + except IndexError as exc: + raise NeedMoreBandsError() from exc else: bands = np.sort(bands) @@ -445,11 +482,8 @@ def nint(num): lumo = [ i[number_electrons / number_electrons_per_band] for i in bands ] # take the n+1th level - except IndexError: - raise ValueError( - "To understand if it is a metal or insulator, " - "need more bands than n_band=number_electrons" - ) + except IndexError as exc: + raise NeedMoreBandsError() from exc if number_electrons % 2 == 1 and len(stored_bands.shape) == 2: # if #electrons is odd and we have a non spin polarized calculation @@ -476,15 +510,9 @@ def nint(num): max_mins = [(max(i), min(i)) for i in levels] if fermi_energy > bands.max(): - raise ValueError( - "The Fermi energy is above all band energies, " - "don't know what to do" - ) + raise FermiEnergyAndBandsEnergiesError(where="above") if fermi_energy < bands.min(): - raise ValueError( - "The Fermi energy is below all band energies, " - "don't know what to do." - ) + raise FermiEnergyAndBandsEnergiesError(where="below") # one band is crossed by the fermi energy if any(i[1] < fermi_energy and fermi_energy < i[0] for i in max_mins): @@ -498,16 +526,14 @@ def nint(num): return False, 0.0, None, None # insulating case else: - # take the max of the band maxima below the fermi energy + # Take the max of the band maxima below the fermi energy. homo = max([i[0] for i in max_mins if i[0] < fermi_energy]) - # take the min of the band minima above the fermi energy + # Take the min of the band minima above the fermi energy.x lumo = min([i[1] for i in max_mins if i[1] > fermi_energy]) gap = lumo - homo if gap <= 0.0: - raise Exception( - "Something wrong has been implemented. " "Revise the code!" - ) + raise WrongCodeError() return True, gap, homo, lumo def search(self, do_all=False): @@ -557,7 +583,7 @@ def search(self, do_all=False): filters["extras.formula"] = {"in": formula_list} if len(self.text_description.value) > 1: - filters["description"] = {"like": f"%{text_description.value}%"} + filters["description"] = {"like": f"%{self.text_description.value}%"} def add_range_filter(bounds, label): filters["extras." + label] = {"and": [{">=": bounds[0]}, {"<": bounds[1]}]} @@ -573,7 +599,7 @@ def add_range_filter(bounds, label): qb.append(WorkChainNode, filters=filters) qb.order_by({WorkChainNode: {"ctime": "desc"}}) - for i, node_tuple in enumerate(qb.iterall()): + for node_tuple in qb.all(): node = node_tuple[0] thumbnail = node.get_extra("thumbnail") node.get_extra("structure_description") diff --git a/nanoribbon/viewers/show_computed.py b/nanoribbon/viewers/show_computed.py index 2ee0d47..d2479f5 100644 --- a/nanoribbon/viewers/show_computed.py +++ b/nanoribbon/viewers/show_computed.py @@ -59,7 +59,7 @@ class BandsViewerWidget(ipw.VBox): selected_band = Int(allow_none=True) selected_kpoint = Int(allow_none=True) selected_spin = Int(allow_none=True) - selected_3D = Int(allow_none=True) + selected_3d = Int(allow_none=True) def __init__(self, **kwargs): self.bands = kwargs["bands"] @@ -111,7 +111,7 @@ def __init__(self, **kwargs): spin_selector = ipw.RadioButtons( options=[("up", 0), ("down", 1)], description="Select spin", disabled=False ) - view_3D = ipw.RadioButtons( + view_3d = ipw.RadioButtons( options=[("no", 0), ("yes", 1)], description="plot3D", disabled=False ) @@ -120,7 +120,7 @@ def __init__(self, **kwargs): band_selector, kpoint_slider, spin_selector, - view_3D, + view_3d, ] plots = [] @@ -134,7 +134,7 @@ def __init__(self, **kwargs): dlink((kpoint_slider, "value"), (self, "selected_kpoint")) dlink((band_selector, "value"), (self, "selected_band")) dlink((spin_selector, "value"), (self, "selected_spin")) - dlink((view_3D, "value"), (self, "selected_3D")) + dlink((view_3d, "value"), (self, "selected_3d")) # Display the orbital map also initially. self.on_band_change(_=None) @@ -245,9 +245,7 @@ def igor_bands(self, ispin): with testio as fobj: fobj.write("IGOR\r") fobj.write("WAVES") - fobj.write( - "\tx1" + ("\ty{}" * nbands).format(*[x for x in range(nbands)]) + "\r" - ) + fobj.write("\tx1" + ("\ty{}" * nbands).format(*range(nbands)) + "\r") fobj.write("BEGIN\r") for i in range(nkpoints): fobj.write(f"\t{k_axis[i]:.7f}") # first column k_axis @@ -737,7 +735,7 @@ def __init__(self, workcalc, **kwargs): self.on_kpoint_change, names=["selected_band", "selected_kpoint", "selected_spin"], ) - self.bands_viewer.observe(self.on_view_3D_change, names="selected_3D") + self.bands_viewer.observe(self.on_view_3d_change, names="selected_3d") # Custom cmap for orbital 2d viewer custom_cmap_colors = [ @@ -774,16 +772,13 @@ def __init__(self, workcalc, **kwargs): self.spin_view_3D.observe(self.on_spin_view_mode_change, names="value") - spin_density_vbox.children += tuple( - [ipw.HTML(value="

Spin density

")] - ) - spin_density_vbox.children += tuple([self.spin_view_3D]) - spin_density_vbox.children += tuple([self.output_s]) + spin_density_vbox.children += (ipw.HTML(value="

Spin density

"),) + spin_density_vbox.children += (self.spin_view_3D,) + spin_density_vbox.children += (self.output_s,) self.on_spin_view_mode_change() - # --------------------------------------- - # STS mapping widget + # STS mapping widget. self.sts_heading = ipw.HTML(value="

LDOS mappings

") self.sts_energy_text = ipw.FloatText( value=round(self._workcalc.get_extra("homo"), 2), description="energy [eV]:" @@ -843,9 +838,9 @@ def _set_viewer_cube_data(viewer): display(self.spinden_viewer_3d) _set_viewer_cube_data(self.spinden_viewer_3d) - def on_view_3D_change(self, _=None): + def on_view_3d_change(self, _=None): """Plot 3D orbitals in case of selection.""" - if self.bands_viewer.selected_3D: + if self.bands_viewer.selected_3d: self.twod_3D = [self.orbital_viewer_2d, self.orbital_viewer_3d] else: self.twod_3D = [self.orbital_viewer_2d] @@ -922,7 +917,7 @@ def on_kpoint_change(self, _=None): with self.output: clear_output() - if self.bands_viewer.selected_3D: + if self.bands_viewer.selected_3d: self.orbital_viewer_3d.arraydata = arraydata hbox = ipw.HBox([self.orbital_viewer_3d]) else: diff --git a/nanoribbon/viewers/utils.py b/nanoribbon/viewers/utils.py index cb66791..d5abbd8 100644 --- a/nanoribbon/viewers/utils.py +++ b/nanoribbon/viewers/utils.py @@ -88,7 +88,7 @@ def from_cube_to_arraydata(cube_content): def adjust_lightness(color, amount=0.9): try: c = mc.cnames[color] - except: + except Exception: c = color c = colorsys.rgb_to_hls(*mc.to_rgb(c)) return colorsys.hls_to_rgb(c[0], max(0, min(1, amount * c[1])), c[2]) From a966e186e21414754746125899ba6c8115930cfc Mon Sep 17 00:00:00 2001 From: Aliaksandr Yakutovich Date: Wed, 31 May 2023 10:59:00 +0000 Subject: [PATCH 3/4] Fix c419 that prevents short-circuiting --- nanoribbon/viewers/search.py | 2 +- nanoribbon/viewers/show_computed.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/nanoribbon/viewers/search.py b/nanoribbon/viewers/search.py index ecc5e2b..0f47e1f 100644 --- a/nanoribbon/viewers/search.py +++ b/nanoribbon/viewers/search.py @@ -222,7 +222,7 @@ def get_calc_by_label(workcalc, label): ] # magnetization ? - if any([k.name[-1].isdigit() for k in structure.kinds]): + if any(k.name[-1].isdigit() for k in structure.kinds): all_steps.append("export_spinden") for label in all_steps: diff --git a/nanoribbon/viewers/show_computed.py b/nanoribbon/viewers/show_computed.py index d2479f5..338e5a5 100644 --- a/nanoribbon/viewers/show_computed.py +++ b/nanoribbon/viewers/show_computed.py @@ -702,7 +702,7 @@ def __init__(self, workcalc, **kwargs): (k, v) for k, v in dict(orbitals_calc.outputs.output_data_multiple).items() ] - elif any(["output_data_multiple" in x for x in orbitals_calc.outputs]): + elif any("output_data_multiple" in x for x in orbitals_calc.outputs): self.list_of_calcs += [ (x, orbitals_calc) for x in orbitals_calc.outputs ] From 34c81e393e206c097f7c81b3b6b372e561555396 Mon Sep 17 00:00:00 2001 From: Aliaksandr Yakutovich Date: Wed, 31 May 2023 12:41:49 +0000 Subject: [PATCH 4/4] Fix remaining pre-commit issue. --- nanoribbon/viewers/cdxml2gnr.py | 38 ++++++++++++++++----------------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/nanoribbon/viewers/cdxml2gnr.py b/nanoribbon/viewers/cdxml2gnr.py index 586fee7..9728a8e 100644 --- a/nanoribbon/viewers/cdxml2gnr.py +++ b/nanoribbon/viewers/cdxml2gnr.py @@ -258,26 +258,24 @@ def _on_file_upload(self, change=None): self.create_cell_btn.disabled = True listmols = [] molid = 0 - for fname, item in change["new"].items(): - frmt = fname.split(".")[-1] - if frmt == "cdxml": - cdxml_file_string = self.file_upload.value[fname]["content"].decode( - "ascii" - ) - self.mols = re.findall( - "") - self.mols[molid] = m - listmols.append( - (str(molid) + ": " + m.formula, molid) - ) # m must be a pb object!!! - molid += 1 - self.allmols.options = listmols - self.allmols.value = 0 - self.allmols.disabled = False - break + + fname = list(change["new"].keys())[0] + frmt = fname.split(".")[-1] + if frmt == "cdxml": + cdxml_file_string = self.file_upload.value[fname]["content"].decode("ascii") + self.mols = re.findall( + "") + self.mols[molid] = m + listmols.append( + (str(molid) + ": " + m.formula, molid) + ) # m must be a pb object!!! + molid += 1 + self.allmols.options = listmols + self.allmols.value = 0 + self.allmols.disabled = False def _on_sketch_selected(self, change=None): self.structure = None # needed to empty view in second viewer