From a6e268ec198c7b3d253e59065a20f17a2db6b87b Mon Sep 17 00:00:00 2001 From: Adam Fekete Date: Mon, 22 Oct 2018 15:06:53 +0100 Subject: [PATCH] adding helper functions to quip tutorial --- quip-gap/scripts/Preprocess.py | 130 +++++++++++++++++++ quip-gap/scripts/Reader.py | 197 +++++++++++++++++++++++++++++ quip-gap/scripts/Visualise.py | 154 ++++++++++++++++++++++ quip-gap/scripts/Visualise_quip.py | 154 ++++++++++++++++++++++ quip-gap/scripts/__init__.py | 0 5 files changed, 635 insertions(+) create mode 100644 quip-gap/scripts/Preprocess.py create mode 100644 quip-gap/scripts/Reader.py create mode 100644 quip-gap/scripts/Visualise.py create mode 100644 quip-gap/scripts/Visualise_quip.py create mode 100644 quip-gap/scripts/__init__.py diff --git a/quip-gap/scripts/Preprocess.py b/quip-gap/scripts/Preprocess.py new file mode 100644 index 0000000..6b29c47 --- /dev/null +++ b/quip-gap/scripts/Preprocess.py @@ -0,0 +1,130 @@ +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() + + diff --git a/quip-gap/scripts/Reader.py b/quip-gap/scripts/Reader.py new file mode 100644 index 0000000..cf9a1ab --- /dev/null +++ b/quip-gap/scripts/Reader.py @@ -0,0 +1,197 @@ +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 diff --git a/quip-gap/scripts/Visualise.py b/quip-gap/scripts/Visualise.py new file mode 100644 index 0000000..1ed65c2 --- /dev/null +++ b/quip-gap/scripts/Visualise.py @@ -0,0 +1,154 @@ +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 diff --git a/quip-gap/scripts/Visualise_quip.py b/quip-gap/scripts/Visualise_quip.py new file mode 100644 index 0000000..1ed65c2 --- /dev/null +++ b/quip-gap/scripts/Visualise_quip.py @@ -0,0 +1,154 @@ +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 diff --git a/quip-gap/scripts/__init__.py b/quip-gap/scripts/__init__.py new file mode 100644 index 0000000..e69de29 -- GitLab