Source code for Cat_structure

# -*- coding: utf-8 -*-
"""
Created on Fri May 26 13:41:45 2017

@author: mpnun
"""

import os
import numpy as np
import copy
import random
import matplotlib.pyplot as plt
import matplotlib as mat
import matplotlib.ticker as mtick

import sys
sys.path.append('C:\Users\mpnun\OneDrive\Documents\\ase')
from ase.build import fcc111, fcc100
from ase.neighborlist import NeighborList
from ase.io import read
from ase.io import write

from metal import metal
from ORR import ORR_rate
from graph_theory import Graph

[docs]class cat_structure(): ''' Catalyst structure with defects ''' def __init__(self, met_name = None, facet = None, dim1 = None, dim2 = None): self.metal = None self.atoms_obj_template = None self.active_atoms = None self.variable_atoms = None self.template_graph = None self.defected_graph = None self.surface_area = None # surface area of the slab in square angstroms self.evaluated = False self.current_density = None self.surf_eng = None self.active_CN = 9 # CN must be less than or equal to this to be active if not met_name is None: ''' Build the slab ''' self.metal = metal(met_name) if facet == '111' or facet == 111: self.atoms_obj_template = fcc111(met_name, size=(dim1, dim2, 4), vacuum=15.0) self.active_CN = 9 elif facet == '100' or facet == 100: self.atoms_obj_template = fcc100(met_name, size=(dim1, dim2, 4), vacuum=15.0) self.active_CN = 8 else: raise ValueError(str(facet) + ' is not a valid facet.') self.atoms_obj_template.set_pbc([True, True, False]) # Find neighbors based on distances rad_list = ( 2.77 + 0.2 ) / 2 * np.ones(len(self.atoms_obj_template)) # list of neighboradii for each site neighb_list = NeighborList(rad_list, self_interaction = False) # set bothways = True to include both ways neighb_list.build(self.atoms_obj_template) self.template_graph = Graph() for i in range(len(self.atoms_obj_template)): self.template_graph.add_vertex(i) for i in range(len(self.atoms_obj_template)): for j in neighb_list.neighbors[i]: self.template_graph.add_edge([i,j]) self.active_atoms = range(2 * dim1 * dim2, 4 * dim1 * dim2) self.variable_atoms = range(3 * dim1 * dim2, 4 * dim1 * dim2) self.defected_graph = copy.deepcopy(self.template_graph) # Compute surface area for use in normalization self.surface_area = np.linalg.norm( np.cross( self.atoms_obj_template.get_cell()[0,:], self.atoms_obj_template.get_cell()[1,:] ) )
[docs] def copy_data(self): ''' Create a new individual with identical essential data ''' child = cat_structure() child.metal = self.metal child.active_atoms = self.active_atoms child.variable_atoms = self.variable_atoms child.template_graph = self.template_graph child.defected_graph = self.defected_graph.copy_data() child.surface_area = self.surface_area child.atoms_obj_template = self.atoms_obj_template return child
[docs] def randomize(self, coverage = None): ''' Randomize the occupancies in the top layer ''' x_rand = [0 for i in range(len(self.variable_atoms)) ] if coverage is None: coverage = random.random() n_occupancies = int( round( coverage * len(self.variable_atoms) ) ) occupied_sites = random.sample(range(len(x_rand)), n_occupancies) for i in occupied_sites: x_rand[i] = 1 return x_rand
[docs] def load_from_file(self, ftoread, d_cut = 0.001): ''' Determine which atoms in the template are missing in the defected structure d_cut: distance cutoff in angstroms ''' ASE_defected = read(ftoread, format = 'xsd') for atom_ind in range(len(self.atoms_obj_template)): # Get position and atomic number of the template atom we are trying to find cart_coords = self.atoms_obj_template.get_positions()[atom_ind, :] atomic_num = self.atoms_obj_template.get_atomic_numbers()[atom_ind] defect_ind = 0 # index of defect atom which might be a match dist = 1.0 # distance between atoms we are trying to match match_found = False while (not match_found) and defect_ind < len(ASE_defected): defect_coords = ASE_defected.get_positions()[defect_ind, :] defect_an = ASE_defected.get_atomic_numbers()[defect_ind] dist = np.linalg.norm( cart_coords - defect_coords ) match_found = (dist < d_cut) #and (defect_an == atomic_num) # We do not need to check whether the elements match defect_ind += 1 if not match_found: self.defected_graph.remove_vertex(atom_ind) self.evaluated = False
[docs] def get_Nnn(self): ''' For each active atom, print the number of nearest neighbors that are also active ''' atom_graph = self.defected_graph for i in self.active_atoms: if atom_graph.is_node(i): if atom_graph.get_coordination_number(i) <= self.active_CN: gcn = atom_graph.get_generalized_coordination_number(i, 12) Nnn = 0 for j in atom_graph.get_neighbors(i): if j in self.active_atoms: if atom_graph.is_node(j): if atom_graph.get_coordination_number(j) <= self.active_CN: Nnn += 1 print [gcn, Nnn]
[docs] def eval_current_density(self, atom_graph = None, normalize = True): ''' Normalized: current density [mA/cm^2] Not normalized: current [mA] ''' if atom_graph is None: atom_graph = self.defected_graph curr = 0 for i in self.active_atoms: if atom_graph.is_node(i): if atom_graph.get_coordination_number(i) <= self.active_CN: gcn = atom_graph.get_generalized_coordination_number(i, 12) BE_OH = self.metal.get_OH_BE(gcn) BE_OOH = self.metal.get_OOH_BE(gcn) curr += ORR_rate(BE_OH, BE_OOH) if normalize: curr = curr / ( self.surface_area * 1.0e-16) # normalize by surface area (in square centimeters) return curr
[docs] def eval_surface_energy(self, atom_graph = None, normalize = True): ''' Normalized: surface energy [J/m^2] Not normalized: formation energy [eV] ''' if atom_graph is None: atom_graph = self.defected_graph E_form = 0 for i in self.active_atoms: if atom_graph.is_node(i): E_form += self.metal.E_coh * ( 1 - np.sqrt( atom_graph.get_coordination_number(i) / 12.0 ) ) #E_form += self.metal.E_coh * ( 1 - atom_graph.get_coordination_number(i) / 12.0 ) if normalize: E_form = E_form * 1.60218e-19 # convert energy from eV to Joules E_form = E_form / ( self.surface_area * 1.0e-20) # normalize by surface area (in square meters) return E_form
[docs] def get_defected_mols(self): ''' Use defected graph and template atoms object to generate the atoms object for the defected structure ''' atoms_obj = copy.deepcopy(self.atoms_obj_template) delete_these = [False for i in range(len(atoms_obj))] for i in self.variable_atoms: if not self.defected_graph.is_node(i): delete_these[i] = True delete_these = np.array(delete_these) del atoms_obj[delete_these] return atoms_obj
[docs] def flip_atom(self, ind): ''' If atom number ind is present in the defected graph, remove it. If it is not present, add it and all edges to adjacent atoms. ''' if self.defected_graph.is_node(ind): self.defected_graph.remove_vertex(ind) else: self.defected_graph.add_vertex(ind) for neighb in self.template_graph.get_neighbors(ind): if self.defected_graph.is_node(neighb): self.defected_graph.add_edge([ind, neighb])
[docs] def show(self, x = None, n_struc = 1, fmat = 'picture'): ''' Print image of surface ''' # Build defected graph if not x is None: self.defected_graph = self.template_graph.copy_data() for i in range(len(x)): if x[i] == 0: self.defected_graph.remove_vertex(self.variable_atoms[i]) # Build ASE atoms object from defected graph defect_atoms_obj = self.get_defected_mols() coords = defect_atoms_obj.get_positions() a_nums = defect_atoms_obj.get_atomic_numbers() chem_symbs = defect_atoms_obj.get_chemical_symbols() # Change top layer atoms to Ni top_layer_z = np.max(coords[:,2]) for atom_ind in range(len(defect_atoms_obj)): if coords[atom_ind,2] > top_layer_z - 0.1: a_nums[atom_ind] = 27 chem_symbs[atom_ind] = 'Co' defect_atoms_obj.set_atomic_numbers(a_nums) defect_atoms_obj.set_chemical_symbols(chem_symbs) if fmat == 'picture': write('structure_' + str(n_struc) + '.png', defect_atoms_obj ) elif fmat == 'xsd': defect_atoms_obj.set_pbc(True) write('structure_' + str(n_struc) + '.xsd', defect_atoms_obj, format = fmat )