adding helper functions to quip tutorial

quip-gap/scripts/Preprocess.py

+from pathlib import Path
+from pprint import pprint
+import json
+from ase.io import read, write
+from ase.geometry import crystal_structure_from_cell
+import numpy as np
+# import numpy.linalg as la
+
+import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
+
+
+class Calculation(object):
+    def __init__(self, *args, **kwargs):
+        self.filepath = kwargs.pop('filepath', None)
+        self.parameters = kwargs
+
+    def get_data(self, index=-1):
+        return read(str(self.filepath))
+
+    @classmethod
+    def from_path(cls, path: Path):
+        with (path / 'gb.json').open() as data_file:
+            gb_data = json.load(data_file)
+
+        with (path / 'subgb.json').open() as data_file:
+            subgb_data = json.load(data_file)
+
+        # print(gb_data['angle'])
+
+        filename = subgb_data['name'] + "_traj.xyz"
+        filepath = (path / filename).resolve()
+        # configuration = read(str((path / filename).resolve()), index=index)
+        # # gb = read(str((path / filename).resolve()), index=-1)
+
+
+        # print('{:=^60}'.format(' '+str(path)+' '))
+        #
+        # print('{:-^40}'.format(' gb.json '))
+        # pprint(gb_data)
+        #
+        # print('{:-^40}'.format(' subgb.json '))
+        # pprint(subgb_data)
+
+        # print('{:-^40}'.format(' initial '))
+        #
+        # force_norm = np.linalg.norm(gb_initial.arrays['force'], axis=1)
+        # force_mean = np.mean(force_norm)
+        # force_std = np.std(force_norm)
+        #
+        # print('Force mean: {:f}, std: {:f}'.format(force_mean, force_std))
+        # pprint(gb_initial.calc.results)
+        #
+        # print('{:-^40}'.format(' final '))
+        #
+        # force_norm = np.linalg.norm(gb_final.arrays['force'], axis=1)
+        # force_mean = np.mean(force_norm)
+        # force_std = np.std(force_norm)
+        #
+        # print('Force mean: {:f}, std: {:f}'.format(force_mean, force_std))
+        # pprint(gb_final.calc.results)
+
+
+        return cls(**{**gb_data, **subgb_data}, filepath=filepath)
+
+
+if __name__ == '__main__':
+
+    # Read grain boundary database
+    dirpath = Path('../GB_alphaFe_001')
+
+    calculations = {
+        'tilt': [Calculation.from_path(calc_dir) for calc_dir in (dirpath / 'tilt').iterdir() if calc_dir.is_dir()],
+        'twist': [Calculation.from_path(calc_dir) for calc_dir in (dirpath / 'twist').iterdir() if calc_dir.is_dir()]
+    }
+
+    # potential energy of the perfect crystal according to a specific potential
+    potential_energy_per_atom = -4.01298214176      # alpha-Fe PotBH
+    eV = 1.6021766208e-19
+    Angstrom = 1.e-10
+
+    angles, energies = [], []
+    for calc in sorted(calculations['tilt'], key=lambda item: item.parameters['angle']):
+
+        # E_gb = calc.parameters.get('E_gb', None)
+        #
+        # if E_gb is None:
+        #     print(calc.filepath)
+        #     print(calc.parameters['converged'])
+        # else:
+
+        # energy = 16.02 / (2 * calc.parameters['A'] ) * \
+        #     (E_gb - potential_energy_per_atom * calc.parameters['n_at'])
+
+        if calc.parameters.get('converged', None):
+            # energy = 16.02 / (2 * calc.parameters['A'] ) * \
+            #     (calc.parameters.get('E_gb') - potential_energy_per_atom * calc.parameters['n_at'])
+            #
+            atoms = calc.get_data()
+            cell = atoms.get_cell()
+            A = cell[0, 0] * cell[1, 1]
+
+            energy = (
+                eV / Angstrom**2 /
+                (2 * A) *
+                (atoms.get_potential_energy() - potential_energy_per_atom * len(atoms))
+            )
+
+            write(calc.filepath.name, atoms)
+
+
+            # print(energy)
+            # print(calc.parameters['converged'])
+            # print(data.get_potential_energy())  # data.get_total_energy() == data.get_potential_energy()
+            # energies.append(calc.parameters['E_gb'] - data.get_total_energy())
+            energies.append(energy)
+            angles.append(calc.parameters['angle'] * 180.0 / np.pi)
+        else:
+            print("not converged: ", calc.filepath)
+
+
+    plt.bar(angles, energies)
+
+    # x_smooth = np.linspace(min(angles), max(angles), 1000, endpoint=True)
+    # f = interp1d(angles, energies, kind='cubic')
+    # plt.plot(x_smooth, f(x_smooth), '-')
+
+    plt.show()
+
+

quip-gap/scripts/Reader.py

+from pathlib import Path
+from pprint import pprint
+import json
+from ase.io import read, write
+from ase.geometry import crystal_structure_from_cell
+import numpy as np
+# import numpy.linalg as la
+
+import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
+
+
+class Calculation(object):
+    def __init__(self, *args, **kwargs):
+        self.filepath = kwargs.pop('filepath', None)
+        self.parameters = kwargs
+
+    def get_data(self, index=-1):
+        return read(str(self.filepath), index=index)
+
+    @classmethod
+    def from_path(cls, path: Path, index=-1):
+        with (path / 'gb.json').open() as data_file:
+            gb_data = json.load(data_file)
+
+        with (path / 'subgb.json').open() as data_file:
+            subgb_data = json.load(data_file)
+
+        # print(gb_data['angle'])
+
+        filename = subgb_data['name'] + "_traj.xyz"
+        filepath = (path / filename).resolve()
+        # configuration = read(str((path / filename).resolve()), index=index)
+        # # gb = read(str((path / filename).resolve()), index=-1)
+
+
+        # print('{:=^60}'.format(' '+str(path)+' '))
+        #
+        # print('{:-^40}'.format(' gb.json '))
+        # pprint(gb_data)
+        #
+        # print('{:-^40}'.format(' subgb.json '))
+        # pprint(subgb_data)
+
+        # print('{:-^40}'.format(' initial '))
+        #
+        # force_norm = np.linalg.norm(gb_initial.arrays['force'], axis=1)
+        # force_mean = np.mean(force_norm)
+        # force_std = np.std(force_norm)
+        #
+        # print('Force mean: {:f}, std: {:f}'.format(force_mean, force_std))
+        # pprint(gb_initial.calc.results)
+        #
+        # print('{:-^40}'.format(' final '))
+        #
+        # force_norm = np.linalg.norm(gb_final.arrays['force'], axis=1)
+        # force_mean = np.mean(force_norm)
+        # force_std = np.std(force_norm)
+        #
+        # print('Force mean: {:f}, std: {:f}'.format(force_mean, force_std))
+        # pprint(gb_final.calc.results)
+
+
+        return cls(**{**gb_data, **subgb_data}, filepath=filepath)
+
+
+if __name__ == '__main__':
+
+    # Read grain boundary database
+    dirpath = Path('../GB_alphaFe_001')
+
+    calculations = {
+        'tilt': [Calculation.from_path(calc_dir) for calc_dir in (dirpath / 'tilt').iterdir() if calc_dir.is_dir()],
+        'twist': [Calculation.from_path(calc_dir) for calc_dir in (dirpath / 'twist').iterdir() if calc_dir.is_dir()]
+    }
+
+    # potential energy of the perfect crystal according to a specific potential
+    potential_energy_per_atom = -4.01298214176      # alpha-Fe PotBH
+    eV = 1.6021766208e-19
+    Angstrom = 1.e-10
+
+    angles, energies = [], []
+    for calc in sorted(calculations['tilt'], key=lambda item: item.parameters['angle']):
+        angles.append(calc.parameters['angle'] * 180.0 / np.pi)
+
+
+        energy = 16.02 / (2 * calc.parameters['A'] ) * \
+            (calc.parameters['E_gb'] - potential_energy_per_atom * calc.parameters['n_at'])
+
+        atoms = calc.get_data()
+        cell = atoms.get_cell()
+        A = cell[0, 0] * cell[1, 1]
+
+        energy = (
+            eV / Angstrom**2 /
+            (2 * A) *
+            (atoms.get_total_energy() - potential_energy_per_atom * len(atoms))
+        )
+
+        print(energy)
+        # print(data.get_potential_energy())  # data.get_total_energy() == data.get_potential_energy()
+        # energies.append(calc.parameters['E_gb'] - data.get_total_energy())
+        energies.append(energy)
+
+
+    plt.bar(angles, energies)
+
+    # x_smooth = np.linspace(min(angles), max(angles), 1000, endpoint=True)
+    # f = interp1d(angles, energies, kind='cubic')
+    # plt.plot(x_smooth, f(x_smooth), '-')
+
+    plt.show()
+
+
+# 00171501160
+# {
+#   "A": 128.39243418897394,
+#   "gbid": "00171501160",
+#   "sigma_csl": 0,
+#   "H": 181.47714791126677,
+#   "n_at": 514,
+#   "orientation_axis": [
+#     0,
+#     0,
+#     1
+#   ],
+#   "boundary_plane": [
+#     1.0,
+#     16.0,
+#     0.0
+#   ],
+#   "coincident_sites": 2,
+#   "angle": 0.12479104151759456,
+#   "type": "symmetric tilt boundary",
+#   "zplanes": [
+#     45.2801,
+#     136.0193
+#   ]
+# }
+# {
+#   "A": 1540.7092103178,
+#   "E_gb": -98740.56229920377,
+#   "gbid": "00171501160_v6bxv2_tv0.4bxv0.1_d1.4z",
+#   "name": "00171501160_v6bxv2_tv0.4bxv0.1_d1.4z",
+#   "area": 1540.7092103178,
+#   "H": 181.47714791,
+#   "n_at": 24636,
+#   "rbt": [
+#     0.4,
+#     0.1
+#   ],
+#   "rcut": 1.4,
+#   "converged": true,
+#   "param_file": "PotBH.xml",
+#   "E_gb_init": -97920.0696568818
+# }
+
+# PotBH
+#
+# Min. Energy: 0.6408 J/m^{2}
+# Min. Energy Structure
+# Microscopic Degrees of Freedom (x,y, Cutoff Criterion)
+# (0.4, 0.1, 1.4)
+# (0.4, 0.1, 1.5)
+# Max. Energy:1.5965 J/m^{2}
+# Max. Energy Structure
+# Microscopic Degrees of Freedom (x,y, Cutoff Criterion)
+# (0.0, 0.1, 2.1)
+
+#
+#
+# models.py
+# gb_ener = 16.02*((sub_dict['E_gb']-(-4.2731*sub_dict['n_at']))/(2*sub_dict['A']))
+#
+# gb_ener = 16.02*((gb_dict['E_gb']-(ener_per_atom[param_file]*float(gb_dict['n_at'])))/(2*gb_dict['A']))
+#
+# calc_gap.py
+# subgbs = [(16.02*(subgb['E_gb']-float(subgb['n_at']*ener_per_atom['PotBH.xml']))/(2.0*subgb['area']), subgb) for subgb in subgbs]
+#
+# calc_energy.py
+# def calc_e_gb(at, E_bulk):
+#   cell = at.get_cell()
+#   A    = cell[0,0]*cell[1,1]
+#   E_gb = (at.get_potential_energy()-(at.n*(E_bulk)))/(2.*A)
+#   print A
+#   print at.get_potential_energy()
+#   print E_gb, 'eV/A^2'
+#   E_gb = 16.02*(at.get_potential_energy()-(at.n*(E_bulk)))/(2.*A)
+    #   print E_gb, 'J/m^2'
+#   return E_gb
+#
+# relax.py
+# E_gb_init   = grain.get_potential_energy()
+# E_gb    = grain.get_potential_energy()
+#
+#
+# ener_per_atom = -4.01298214176
\ No newline at end of file

quip-gap/scripts/Visualise.py

+import io
+import uuid
+
+from nglview import register_backend, Structure
+from ipywidgets import Dropdown, FloatSlider, IntSlider, HBox, VBox, Output
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+    
+
+@register_backend('ase')
+class MyASEStructure(Structure):
+    def __init__(self, atoms, bfactor=[], occupancy=[]):
+        # super(MyASEStructure, self).__init__()
+        self.ext = 'pdb'
+        self.params = {}
+        self._atoms = atoms
+        self.bfactor = bfactor  # [min, max]
+        self.occupancy = occupancy  # [0, 1]
+        self.id = str(uuid.uuid4())
+
+    def get_structure_string(self):
+        """ PDB file format:
+        CRYST1   16.980   62.517  124.864  90.00  90.00  90.00 P 1
+        MODEL     1
+        ATOM      0   Fe MOL     1      15.431  60.277   6.801  1.00  0.00          FE
+        ATOM      1   Fe MOL     1       1.273   3.392  93.940  1.00  0.00          FE
+        """
+        # with io.StringIO() as stream:
+        #     self.structure.write(stream, format='proteindatabank')
+        #     data = stream.getvalue()
+
+        data = ""
+
+        if self._atoms.get_pbc().any():
+            cellpar = self._atoms.get_cell_lengths_and_angles()
+
+            str_format = 'CRYST1' + '{:9.3f}' * 3 + '{:7.2f}' * 3 + ' P 1\n'
+            data += str_format.format(*cellpar.tolist())
+
+        data += 'MODEL     1\n'
+
+        str_format = 'ATOM  {:5d} {:>4s} MOL     1    {:8.3f}{:8.3f}{:8.3f}{:6.2f}{:6.2f}          {:2s}\n'
+        for index, atom in enumerate(self._atoms):
+            data += str_format.format(
+                index,
+                atom.symbol,
+                atom.position[0].tolist(),
+                atom.position[1].tolist(),
+                atom.position[2].tolist(),
+                self.occupancy[index] if index <= len(self.occupancy) - 1 else 1.0,
+                self.bfactor[index] if index <= len(self.bfactor) - 1 else 1.0,
+                atom.symbol.upper()
+            )
+
+        data += 'ENDMDL\n'
+
+        return data
+
+
+def ViewStructure(atoms):
+    import nglview
+
+    view = nglview.NGLWidget()
+    
+    structure = MyASEStructure(atoms)
+    view.add_structure(structure)
+    
+    return view
+
+
+
+class AtomViewer(object):
+    def __init__(self, atoms, data=[], xsize=1000, ysize=500):
+        self.view = self._init_nglview(atoms, data, xsize, ysize)
+
+        self.widgets = {
+            'radius': FloatSlider(
+                value=0.8, min=0.0, max=1.5, step=0.01,
+                description='Ball size'
+            ),
+            'color_scheme': Dropdown(description='Solor scheme:'),
+            'colorbar': Output()
+        }
+        self.show_colorbar(data)
+
+        self.widgets['radius'].observe(self._update_repr)
+
+        self.gui = VBox([
+            self.view,
+            self.widgets['colorbar'],
+            self.widgets['radius']])
+
+    def _update_repr(self, chg=None):
+        self.view.update_spacefill(
+            radiusType='radius',
+            radius=self.widgets['radius'].value
+        )
+
+    def show_colorbar(self, data):
+        with self.widgets['colorbar']:
+            # Have colormaps separated into categories:
+            # http://matplotlib.org/examples/color/colormaps_reference.html
+            cmap = 'rainbow'
+
+            fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 2))
+            img = ax1.imshow([[min(data), max(data)]], aspect='auto', cmap=plt.get_cmap(cmap))
+            ax1.remove()
+            cbar = fig.colorbar(img, cax=ax2, orientation='horizontal')
+
+            plt.show()
+
+    @staticmethod
+    def _init_nglview(atoms, data, xsize, ysize):
+        import nglview
+
+        view = nglview.NGLWidget(gui=False)
+        view._remote_call(
+            'setSize',
+            target='Widget',
+            args=[
+                '{:d}px'.format(xsize),
+                '{:d}px'.format(ysize)
+            ]
+        )
+
+        data = np.max(data)-data
+
+        structure = MyASEStructure(atoms, bfactor=data)
+        view.add_structure(structure)
+
+        view.clear_representations()
+        view.add_unitcell()
+
+        view.add_spacefill(
+            # radiusType='radius',
+            # radius=1.0,
+            color_scheme='bfactor',
+            color_scale='rainbow'
+        )
+        view.update_spacefill(
+            radiusType='radius',
+            radius=1.0
+        )
+
+
+        # update camera type
+        view.control.spin([1, 0, 0], np.pi / 2)
+        view.control.spin([0, 0, 1], np.pi / 2)
+        view.camera = 'orthographic'
+        view.center()
+
+        return view

quip-gap/scripts/Visualise_quip.py

+import io
+import uuid
+
+from nglview import register_backend, Structure
+from ipywidgets import Dropdown, FloatSlider, IntSlider, HBox, VBox, Output
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+    
+
+@register_backend('ase')
+class MyASEStructure(Structure):
+    def __init__(self, atoms, bfactor=[], occupancy=[]):
+        # super(MyASEStructure, self).__init__()
+        self.ext = 'pdb'
+        self.params = {}
+        self._atoms = atoms
+        self.bfactor = bfactor  # [min, max]
+        self.occupancy = occupancy  # [0, 1]
+        self.id = str(uuid.uuid4())
+
+    def get_structure_string(self):
+        """ PDB file format:
+        CRYST1   16.980   62.517  124.864  90.00  90.00  90.00 P 1
+        MODEL     1
+        ATOM      0   Fe MOL     1      15.431  60.277   6.801  1.00  0.00          FE
+        ATOM      1   Fe MOL     1       1.273   3.392  93.940  1.00  0.00          FE
+        """
+        # with io.StringIO() as stream:
+        #     self.structure.write(stream, format='proteindatabank')
+        #     data = stream.getvalue()
+
+        data = ""
+
+        if self._atoms.get_pbc().any():
+            cellpar = self._atoms.get_cell_lengths_and_angles()
+
+            str_format = 'CRYST1' + '{:9.3f}' * 3 + '{:7.2f}' * 3 + ' P 1\n'
+            data += str_format.format(*cellpar.tolist())
+
+        data += 'MODEL     1\n'
+
+        str_format = 'ATOM  {:5d} {:>4s} MOL     1    {:8.3f}{:8.3f}{:8.3f}{:6.2f}{:6.2f}          {:2s}\n'
+        for index, atom in enumerate(self._atoms):
+            data += str_format.format(
+                index,
+                atom.symbol,
+                atom.position[0].tolist(),
+                atom.position[1].tolist(),
+                atom.position[2].tolist(),
+                self.occupancy[index] if index <= len(self.occupancy) - 1 else 1.0,
+                self.bfactor[index] if index <= len(self.bfactor) - 1 else 1.0,
+                atom.symbol.upper()
+            )
+
+        data += 'ENDMDL\n'
+
+        return data
+
+
+def ViewStructure(atoms):
+    import nglview
+
+    view = nglview.NGLWidget()
+    
+    structure = MyASEStructure(atoms)
+    view.add_structure(structure)
+    
+    return view
+
+
+
+class AtomViewer(object):
+    def __init__(self, atoms, data=[], xsize=1000, ysize=500):
+        self.view = self._init_nglview(atoms, data, xsize, ysize)
+
+        self.widgets = {
+            'radius': FloatSlider(
+                value=0.8, min=0.0, max=1.5, step=0.01,
+                description='Ball size'
+            ),
+            'color_scheme': Dropdown(description='Solor scheme:'),
+            'colorbar': Output()
+        }
+        self.show_colorbar(data)
+
+        self.widgets['radius'].observe(self._update_repr)
+
+        self.gui = VBox([
+            self.view,
+            self.widgets['colorbar'],
+            self.widgets['radius']])
+
+    def _update_repr(self, chg=None):
+        self.view.update_spacefill(
+            radiusType='radius',
+            radius=self.widgets['radius'].value
+        )
+
+    def show_colorbar(self, data):
+        with self.widgets['colorbar']:
+            # Have colormaps separated into categories:
+            # http://matplotlib.org/examples/color/colormaps_reference.html
+            cmap = 'rainbow'
+
+            fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 2))
+            img = ax1.imshow([[min(data), max(data)]], aspect='auto', cmap=plt.get_cmap(cmap))
+            ax1.remove()
+            cbar = fig.colorbar(img, cax=ax2, orientation='horizontal')
+
+            plt.show()
+
+    @staticmethod
+    def _init_nglview(atoms, data, xsize, ysize):
+        import nglview
+
+        view = nglview.NGLWidget(gui=False)
+        view._remote_call(
+            'setSize',
+            target='Widget',
+            args=[
+                '{:d}px'.format(xsize),
+                '{:d}px'.format(ysize)
+            ]
+        )
+
+        data = np.max(data)-data
+
+        structure = MyASEStructure(atoms, bfactor=data)
+        view.add_structure(structure)
+
+        view.clear_representations()
+        view.add_unitcell()
+
+        view.add_spacefill(
+            # radiusType='radius',
+            # radius=1.0,
+            color_scheme='bfactor',
+            color_scale='rainbow'
+        )
+        view.update_spacefill(
+            radiusType='radius',
+            radius=1.0
+        )
+
+
+        # update camera type
+        view.control.spin([1, 0, 0], np.pi / 2)
+        view.control.spin([0, 0, 1], np.pi / 2)
+        view.camera = 'orthographic'
+        view.center()
+
+        return view