# -*- coding: utf-8 -*-
"""
Created on Tue May 19 22:28:09 2020
@author: lansf
"""
from __future__ import absolute_import, division, print_function
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
from ase.io import read
from .vasp_dos import get_band_center
from .coordination import Coordination
from .plotting_tools import set_figure_settings
[docs]def get_adsorbate_indices(GAS_CONTCAR, ADSORBATE_CONTCAR):
""" Identify the adsorbate and adsorption site indices
Parameters
----------
GAS_CONTCAR : str
Location to VASP CONTCAR file for a gas molecule.
ADSORBATE_CONTCAR : str
Location to VASP CONTCAR file for an adsorbate + surface
Returns
-------
adsorbate_indices : list[int]
Index (indices) of adsorbate atom(s) in REFERENCE_PDOS. Obtained
from pdos_overlap.get_adsorbate_indices.
site_indices : list[int]
Index (indices) of adsorbation site atom(s) in REFERENCE_PDOS.
Obtained from pdos_overlap.get_adsorbate_indices.
"""
gas_symbols = list(set(read(GAS_CONTCAR).get_chemical_symbols()))
adsorbate_symbols = read(ADSORBATE_CONTCAR).get_chemical_symbols()
adsorbate_indices = np.arange(len(adsorbate_symbols))[np.isin(adsorbate_symbols,gas_symbols)]
CN = Coordination(read(ADSORBATE_CONTCAR),exclude=adsorbate_indices,cutoff=1.25)
bonded = CN.get_bonded()
site_indices = []
for i in adsorbate_indices:
site_indices += bonded[i]
site_indices = list(set(site_indices))
return adsorbate_indices, site_indices
[docs]class PDOS_OVERLAP:
"""Class for calculating adsorbate-surface relative energy overlap"""
def __init__(self, GAS_PDOS, REFERENCE_PDOS, adsorbate_indices, site_indices\
, min_occupation=0.9, upshift=4.5, energy_weight=3\
, sum_density=False, sum_spin=True):
"""
Parameters
----------
GAS_PDOS : vasp_dos.VASP_DOS
VASP_DOS object of gas phase calculation of adsorbate.
REFERENCE_PDOS : vasp_dos.VASP_DOS
VASP_DOS object of reference adsorbate+surface system.
adsorbate_indices : list[int]
Index (indices) of adsorbate atom(s) in REFERENCE_PDOS. Obtained
from pdos_overlap.get_adsorbate_indices.
site_indices : list[int]
Index (indices) of adsorbation site atom(s) in REFERENCE_PDOS.
Obtained from pdos_overlap.get_adsorbate_indices.
min_occupation : float
Minimum state occupation for a section of density to be considered
a molecular orbital.
upshift : float
The energy the gas molecular orbitals are shifted before pairing
energy_weight : float
The degree to weight the energy error before computing the
associated liklihood that goes into computing the scores. Increasing
the energy weight sharpens the univariate distribution around
matching peak positions.
sum_density : bool
Tag indicates if state density on the same sublevel should be
summed when calculating overlap of the site atomic orbitals and the
adsorbate molecular orbitals.
sum_spin : bool
Tag indictes if state density with different spins should be summed
when generating gas and adsorbate orbital features.
Attributes
----------
gas_indices : list[ing]
Atom index (indices) of the gas molecule.
gas_features : numpy.ndarray
gas molecular orbital features used in matching gas and adsorbate
molecular orbitals
adsorbate_features : numpy.ndarray
Adsorbate molecular orbital features used in matching gas and
adsorbate molecular orbitals
gas_orbitals : numpy.ndarray
Length mxn array where m is the number of molecular orbitals and
n is equal to GAS_PDOS.ndos (dos discritizations)
gas_occupations : numpy.ndarray
Integrated state density of each gas orbital.
adsorbate_orbitals : numpy.ndarray
Length mxn array where m is the number of adsorbate orbitals and
n is equal to REFERENCE_PDOS.ndos (dos discritizations)
adsorbate_occupations : numpy.ndarray
Integrated state density of each adsorbate molecular orbital.
orbital_scores : numpy.ndarray
Array of size mxn where m and n are the number of gas and adsorbate
molecular orbitals respecitvely. Indicates similarity. If
the number of gas orbitals equals the number of adsorbate orbitals
then the scores are orbital-pairing scores.
feature_type : str
Indicates whether integrated state density of atomic orbitals are
used for molecular orbital features or if molecular orbital
occupation fractions.
gas_orbital_indices : list
List of vector indices that aportion the gas molecule energy and
projected density into molecular orbitals
adsorbate_orbital_indices : list
List of vector indices that aportion the adsorbate molecule energy
and projected density into molecular orbitals
gas_band_centers : numpy.ndarray
Gas molecular orbital band centers
adsorbate_band_centers : numpy.ndarray
Adsorbate molecular orbital band centers
gas_2_adsorbate : numpy.ndarray
Array of gas and adsorbate orbital indices along with corresponding
band centers
Notes
-----
GAS_PDOS is used to determine the number of molecular orbitals that can
interact with the surface and to calculate relative pdos overlap the
with projected density of adsorption sites without adsorbates.
REFERENCE_PDOS is used to determine which metal states, projected onto
atomic orbitals, will interact with the adsorbate/gas molecular orbitals.
All parameters are saved as attributes.
"""
self.GAS_PDOS = GAS_PDOS
self.REFERENCE_PDOS = REFERENCE_PDOS
self.adsorbate_indices = adsorbate_indices
self.site_indices = site_indices
self.min_occupation = min_occupation
self.upshift = upshift
self.energy_weight = energy_weight
self.sum_density = sum_density
self.sum_spin = sum_spin
gas_indices = list(range(GAS_PDOS.natoms))
gas_peak_energies, gas_peak_densities, gas_orbitals\
, gas_orbital_indices, gas_occupations\
= self._get_mol_orbitals(GAS_PDOS, gas_indices)
adsorbate_peak_energies, adsorbate_peak_densities, adsorbate_orbitals \
, adsorbate_orbital_indices, adsorbate_occupations \
= self._get_mol_orbitals(REFERENCE_PDOS, adsorbate_indices)
self.adsorbate_orbitals = adsorbate_orbitals
if gas_orbitals.shape[0] == adsorbate_orbitals.shape[0]:
self.feature_type = 'states'
else:
self.feature_type = 'state fractions'
gas_features\
= self._get_orbital_features(GAS_PDOS, gas_orbital_indices\
, gas_indices, upshift=self.upshift)
adsorbate_features\
= self._get_orbital_features(REFERENCE_PDOS, adsorbate_orbital_indices\
, adsorbate_indices, upshift=0)
orbital_scores = self.get_orbital_scores(gas_features\
, adsorbate_features, energy_weight)
self.gas_band_centers = get_band_center(GAS_PDOS.get_energies()\
, gas_orbitals)
self.adsorbate_band_centers = get_band_center(\
REFERENCE_PDOS.get_energies()\
, adsorbate_orbitals)
self.gas_2_adsorbate = self._assign_orbitals(orbital_scores)
self.gas_indices = gas_indices
self.gas_features = gas_features
self.adsorbate_features = adsorbate_features
self.gas_orbitals = gas_orbitals
self.gas_occupations = gas_occupations
self.adsorbate_occupations = adsorbate_occupations
self.orbital_scores = orbital_scores
self.gas_orbital_indices = gas_orbital_indices
self.adsorbate_orbital_indices = adsorbate_orbital_indices
[docs] def _get_mol_orbitals(self, PDOS, atom_indices):
""" Molecular orbitals as an M x ndos array
Parameters
----------
PDOS : vasp_dos.VASP_DOS
VASP_DOS object of gas or surface phase calculation of adsorbate
atom_indices : list[int]
indices of atoms to include in the TOTAL DOS. If empty, then all
atoms are included (as for the gas)
Returns
-------
peak_energies : float or numpy.ndarray
peak energy(s) of molecular orbital(s)
peak_densities : float or numpy.ndarray
peak density(ies) of molecular orbital(s)
mol_orbitals : numpy.ndarray
M x ndos 2D array where M is the number of molecular orbitals and
ndos is the number of gridpoints for the density of states
orbital_indices : list[numpy.ndarray]
Length M list of non-zero orbital indices
"""
min_occupation = self.min_occupation
TOTAL_DOS = np.zeros(PDOS.ndos)
#sum over projected density of states for each atom
orbital_list = [key for key in PDOS.orbital_dictionary.keys()]
orbital_list, projected_dos = PDOS.get_site_dos(atom_indices, orbital_list)
TOTAL_DOS = projected_dos.sum(axis=0)
#obtain all local maxima (including endpoints) - True or False
peaks = np.r_[True, TOTAL_DOS[1:] > TOTAL_DOS[:-1]]\
& np.r_[TOTAL_DOS[:-1] > TOTAL_DOS[1:], True]
#endpoints should be zero density
all_indices = np.arange(PDOS.ndos)
peak_indices = all_indices[peaks]
num_peaks = peak_indices.size
peak_energies = PDOS.get_energies()[peak_indices]
peak_densities = TOTAL_DOS[peak_indices]
if num_peaks == 1:
mol_orbitals = TOTAL_DOS
orbital_indices = all_indices
else:
mol_orbitals = np.zeros((num_peaks,TOTAL_DOS.size))
orbital_indices = []
midpoints = (0.5 * peak_indices[0:-1] + 0.5 * peak_indices[1:]).astype(int)
mol_orbitals[0][0:midpoints[0]] = TOTAL_DOS[0:midpoints[0]]
orbital_indices.append(all_indices[0:midpoints[0]])
for i in range(midpoints.size-1):
mol_orbitals[i+1][midpoints[i]:midpoints[i+1]] = TOTAL_DOS[midpoints[i]:midpoints[i+1]]
orbital_indices.append(all_indices[midpoints[i]:midpoints[i+1]])
mol_orbitals[-1][midpoints[-1]:] = TOTAL_DOS[midpoints[-1]:]
orbital_indices.append(all_indices[midpoints[-1]:])
orbital_occupations = np.trapz(mol_orbitals, PDOS.get_energies(), axis=1)
new_mol_orbitals = []
new_orbital_indices = []
new_occupations = []
orbital_num = 0
#theoreticly new_occupations < 1 but account for some loss of electrons
while orbital_num < num_peaks:
new_mol_orbitals.append(mol_orbitals[orbital_num])
new_orbital_indices.append(orbital_indices[orbital_num])
new_occupations.append(orbital_occupations[orbital_num])
orbital_num += 1
while orbital_num < num_peaks\
and orbital_occupations[orbital_num - 1] < min_occupation\
and orbital_occupations[orbital_num] < min_occupation:
new_mol_orbitals[-1] += mol_orbitals[orbital_num]
new_orbital_indices[-1]\
= np.concatenate( ( new_orbital_indices[-1]\
, orbital_indices[orbital_num] ) )
new_occupations[-1] += orbital_occupations[orbital_num]
orbital_num += 1
num_new_orbitals = len(new_mol_orbitals)
band_centers = get_band_center(PDOS.get_energies(),np.array(new_mol_orbitals))
num_orbitals_final = len([value for value in new_occupations if value > min_occupation])
mol_orbitals_final = np.zeros((num_orbitals_final,TOTAL_DOS.size))
orbital_indices_final = [np.array([]).astype(int) for _ in range(num_orbitals_final)]
occupations_final = np.zeros(num_orbitals_final)
j = 0
for i in range(num_new_orbitals):
if new_occupations[i] > min_occupation or i == 0:
mol_orbitals_final[j] += new_mol_orbitals[i]
orbital_indices_final[j] = np.concatenate(
(orbital_indices_final[j],new_orbital_indices[i]) )
occupations_final[j] += new_occupations[i]
#only update if a complete orbital is found
if new_occupations[i] > min_occupation:
j += 1
elif i == num_new_orbitals - 1:
mol_orbitals_final[-1] += new_mol_orbitals[i]
orbital_indices_final[-1] = np.concatenate(
(orbital_indices_final[-1],new_orbital_indices[i]) )
occupations_final[-1] += new_occupations[i]
elif abs(band_centers[i] - band_centers[i-1])\
<= abs(band_centers[i] - band_centers[i+1]):
mol_orbitals_final[j-1] += new_mol_orbitals[i]
orbital_indices_final[j-1] = np.concatenate(
(orbital_indices_final[j-1],new_orbital_indices[i]) )
occupations_final[j-1] += new_occupations[i]
elif abs(band_centers[i] - band_centers[i-1])\
> abs(band_centers[i] - band_centers[i+1]):
mol_orbitals_final[j] += new_mol_orbitals[i]
orbital_indices_final[j] = np.concatenate(
(orbital_indices_final[j],new_orbital_indices[i]) )
occupations_final[j] += new_occupations[i]
molecular_orbitals = mol_orbitals_final
orbital_indices = orbital_indices_final
occupations = occupations_final
return peak_energies, peak_densities, molecular_orbitals\
, orbital_indices, occupations
[docs] def _get_orbital_features(self, PDOS, orbital_indices, atom_indices, upshift):
""" Molecular orbitals as an M x ndos array
Parameters
----------
PDOS : vasp_dos.VASP_DOS
VASP_DOS object of gas or surface phase calculation of adsorbate
orbital_indices : list[numpy.ndarray]
Length M list of non-zero orbital indices
atom_indices : list[int]
indices of atoms to include in the TOTAL PDOS.
upshift : float
The energy the gas molecular orbitals are shifted before pairing
Returns
-------
orbital_features : numpy.ndarray
M x n_features 2D array where M is the number of molecular orbitals
and n_features is the number of orbital features for calculating
orbital similarity. Includes molecular orbital energy center and
integrated s, p, and d molecular orbital density of states
"""
feature_type = self.feature_type
sum_spin = self.sum_spin
num_orbitals = len(orbital_indices)
energies = PDOS.get_energies() + upshift
orbital_list = list(PDOS.orbital_dictionary.keys())
#feature_list = ['s','py', 'pz', 'px', 'dxy', 'dyz', 'dz2', 'dxz', 'dx2-y2']
feature_list = ['s','py', 'pz', 'px']
if sum_spin == False:
orbital_list=['s+', 's-', 'py+', 'py-', 'pz+', 'pz-', 'px+', 'px-']
num_features = len(feature_list) + 1
TOTAL_DOS = np.zeros(PDOS.ndos)
TOTAL_PDOS = np.zeros((num_features - 1,PDOS.ndos))
#sum over projected density of states for each atom
orbital_list, projected_dos = PDOS.get_site_dos(atom_indices\
, orbital_list, sum_density=False, sum_spin=sum_spin)
TOTAL_DOS = projected_dos.sum(axis=0)
for count, value in enumerate(feature_list):
if value in orbital_list:
TOTAL_PDOS[count] += projected_dos[orbital_list.index(value)]
orbital_features = np.zeros((num_orbitals,num_features))
for count, index_values in enumerate(orbital_indices):
orbital_features[count][0] = get_band_center(energies[index_values],TOTAL_DOS[index_values])
for count2, feature_level in enumerate(TOTAL_PDOS):
if feature_type == 'states':
orbital_features[count][count2+1]\
= np.trapz(feature_level[index_values], energies[index_values])
elif feature_type == 'state fractions':
orbital_features[count][count2+1]\
= np.trapz(feature_level[index_values], energies[index_values])\
/ np.trapz(TOTAL_DOS[index_values], energies[index_values])
return orbital_features
[docs] @staticmethod
def _get_orbital_scores(reference_features, comparison_features\
, energy_weight=3):
""" Orbital scores given reference and comparision features
Parameters
----------
reference_features : numpy.ndarray
Reference features of molecular orbitals whose values are treated
as the mean
comparison_features : numpy.ndarray
Comparison features that are treated as perturbations from the
reference features for generating probabilities
energy_weight : float
The degree to weight the energy error before computing the
associated liklihood that goes into computing the scores. Increasing
the energy weight sharpens the univariate distribution around
matching peak positions.
Returns
-------
orbital_scores : numpy.ndarray
M x N 2D array where M is the number of reference orbitals
and N is the number of comparison orbitals
"""
num_reference_orbitals = reference_features.shape[0]
num_comparison_orbitals = comparison_features.shape[0]
orbital_scores = np.zeros((num_reference_orbitals, num_comparison_orbitals))
for count in range(num_reference_orbitals):
diff_squared = (comparison_features - reference_features[count])**2
var_feature = np.sum(diff_squared,axis=0) / num_comparison_orbitals
var_feature[var_feature[...] <= 10**-6] = 10 #prevent divide by zero error
#gaussian based likelihood
diff_squared[:,0] *= energy_weight
univariate = np.exp(-1 * diff_squared / var_feature)
orbital_scores[count] = np.prod(univariate, axis=1)
return orbital_scores
[docs] @staticmethod
def get_pair_scores(orbital_scores, max_iterations = 500, verbose=False):
""" Get normalized pair scores using matrix scaling
Parameters
----------
orbital_scores : numpy.ndarray
M x M 2D array
max_iterations : int
Maximum number of iterations before ending normalization procedure
verbose : bool
Print completion information
Returns
-------
orbital_scores : numpy.ndarray
M x M 2D array normalized symmetric matrix.
Notes
-----
Matrix scaling perfomed using RAS method described in
DOI: 10.1109/FOCS.2017.87
"""
error = 1
iteration = 0
orbital_scores = orbital_scores.copy()
while error > 10**-7 and iteration < max_iterations:
orbital_scores = (orbital_scores * orbital_scores.T)**0.5
orbital_scores /= orbital_scores.sum(axis=0, keepdims=True)
orbital_scores /= orbital_scores.sum(axis=1, keepdims=True)
error = (abs(1 - orbital_scores.sum(axis=0))).max()
iteration += 1
if verbose == True:
print('The max error is ' + str(error))
print('Number of iterations is ' + str(iteration))
return orbital_scores
[docs] def get_energy_overlap(self, atomic_orbitals
,sum_overlap=False, sum_spin=True, normalize=False):
""" Calculate energy overlap of molecular orbitals with atomic orbitals
Parameters
----------
atomic_orbitals : list[str]
Adsorption site atomic orbitals of interest.
sum_density : bool
Tag indicates if overlap of the site atomic orbitals and the
adsorbate molecular orbitals shold be summed for overlap of atomic
orbitals in the same sublevel.
sum_spin : bool
Tag indictes if state density with different spins should be summed
when generating gas and adsorbate orbital features.
normalize : bool
Tag indicates if energy overlap should be normalized.
Returns
-------
overlap_orbitals : list[str]
List of atomic orbitals for which overlap energies of the adsorbate
state with the metal adsorption-site states are calculated.
energy_overlap : numpy.ndarray
Array of shape nxp where n is the number of adsorbate molecular
orbitals and p is the number of adsorption-site atomic orbitals.
Elements are the projected density of states overlap integral of
an adsorption site states (projected onto atomic orbitals), and
the adsorbate states (projected and summed into molecular orbitals)
Notes
-----
As implemented, molecular_orbitals are for the adsorbate and PDOS is
for some surface.
"""
site_indices = self.site_indices
REFERENCE_PDOS = self.REFERENCE_PDOS
e_fermi = REFERENCE_PDOS.e_fermi
adsorbate_orbitals = np.atleast_2d(self.adsorbate_orbitals.copy())
orbital_list, TOTAL_PDOS = REFERENCE_PDOS.get_site_dos(site_indices\
, atomic_orbitals, sum_density=False, sum_spin=sum_spin)
TOTAL_PDOS = np.atleast_2d(TOTAL_PDOS)
energies = REFERENCE_PDOS.get_energies()
#ensure only populated shates are used
for i in range(adsorbate_orbitals.shape[0]):
adsorbate_orbitals[i][energies[...] > e_fermi] = 0
for i in range(TOTAL_PDOS.shape[0]):
TOTAL_PDOS[i][energies[...] > e_fermi] = 0
#calculate sumj(Mij*projected_dosj) for all i
overlap_orbitals = orbital_list
energy_overlap = np.zeros((adsorbate_orbitals.shape[0], TOTAL_PDOS.shape[0]))
for i in range(adsorbate_orbitals.shape[0]):
if normalize == False:
energy_overlap[i] = np.trapz(adsorbate_orbitals[i]**0.5 * TOTAL_PDOS**0.5\
, energies)
else:
numerator = np.trapz(adsorbate_orbitals[i]**0.5
* TOTAL_PDOS**0.5, energies)
denominator = (np.trapz(adsorbate_orbitals[i], energies)**0.5
* np.trapz(TOTAL_PDOS, energies)**0.5 )
denominator[denominator[...] == 0] = 1 # avoid devide by zero
energy_overlap[i] = numerator / denominator
if sum_overlap == True:
new_overlap_orbitals = []
new_energy_overlap = []
for count, orbital in enumerate(overlap_orbitals):
if orbital[0] not in new_overlap_orbitals:
new_overlap_orbitals.append(orbital[0])
new_energy_overlap.append(energy_overlap[count])
else:
new_energy_overlap[-1] += energy_overlap[count]
overlap_orbitals = np.array(new_overlap_orbitals)
energy_overlap = np.array(new_energy_overlap)
return overlap_orbitals, energy_overlap.squeeze()
[docs] def get_orbital_interaction(self,orbital_index, PDOS, site_indices\
, atomic_orbitals, BULK_PDOS, bulk_atom=0\
, sum_interaction=False, sum_spin=True\
, method='orbital_bond_energy'\
, use_orbital_proximity=False\
, index_type='gas'):
""" Calculate surface and gas orbital interactions
Parameters
----------
orbital_index : int
Identifies the gas/adsorbate orbital with which to calculate the
overlap of the PDOS atomic orbitals
PDOS : vasp_dos.VASP_DOS
VASP_DOS object of surface phase calculation of adsorbate
overlap_orbitals : list[str]
List of atomic orbitals for which overlap energies of the gas state
with a metal adsorption-site states are calculated.
site_indices : list[int]
indices of atoms to include in the TOTAL DOS.
atomic_orbitals : list[str]
Adsorption site atomic orbitals of interest.
sum_interaction : bool
Tag indicates if interactions of molecular orbitals with atomic
orbitals on the same sublevel should be summed
sum_spin : bool
Tag indictes if state density with different spins should be summed
when generating gas and adsorbate orbital features.
method : str
Can be 'orbital_bond_energy' or 'band_width' and dictates whate
values from the pdos are used in computing relative potential to bond
use_orbital_proximity :bool
Indicates whether proximity of the surface to the gas states should
be used to scale interactions. Only used if index_type = 'gas'
index_type : str
'gas' or 'adsorbate' identifies whether the orbital indices are for
the gs or adsorbate orbitals
Returns
-------
orbital_interaction : numpy.ndarray
Metal atomic orbital interactions with the gas molecular orbitals
"""
get_energy_overlap = self.get_energy_overlap
if index_type == 'gas':
#get positions of the gas indices in the conversion matrix
gas_positions = [i for i in range(self.gas_2_adsorbate.shape[0])\
if self.gas_2_adsorbate[i][0] == orbital_index]
gas_energy = self.gas_band_centers[orbital_index]
#get the adsorbate indices
adsorbate_indices = self.gas_2_adsorbate[gas_positions,1].astype('int')
elif index_type == 'adsorbate':
adsorbate_indices = orbital_index
overlap_orbitals, normalized_overlap = get_energy_overlap(atomic_orbitals
, sum_overlap=False
, sum_spin=sum_spin
, normalize=True)
# This sums overlap if index_type == 'gas' and a single gas orbital
# matches more than one adsorabte orbital
energy_overlap = np.array([normalized_overlap[adsorbate_indices\
,overlap_orbitals.index(i)].sum()\
for i in overlap_orbitals])
if use_orbital_proximity == True and index_type == 'gas':
orbital_proximity = PDOS.get_orbital_proximity(gas_energy\
, site_indices, atomic_orbitals\
, sum_density=False, sum_spin=sum_spin)
bulk_proximity = BULK_PDOS.get_orbital_proximity(gas_energy\
, bulk_atom, atomic_orbitals\
, sum_density=False, sum_spin=sum_spin)
scaling_factor = bulk_proximity / orbital_proximity
else:
scaling_factor = 1
if method == 'orbital_bond_energy':
bulk_bond_energy = BULK_PDOS.get_bond_energy(bulk_atom, atomic_orbitals\
, sum_density=False, sum_spin=sum_spin)
bond_energy = PDOS.get_bond_energy(site_indices, atomic_orbitals\
, sum_density=False, sum_spin=sum_spin)
orbital_interaction = (bulk_bond_energy - bond_energy)\
* energy_overlap * scaling_factor
elif method == 'band_width':
bulk_moment = BULK_PDOS.get_second_moment(bulk_atom, atomic_orbitals\
, sum_density=False, sum_spin=sum_spin)
moment = PDOS.get_second_moment(site_indices, atomic_orbitals\
, sum_density=False, sum_spin=sum_spin)
orbital_interaction = (moment**0.5 - bulk_moment**0.5)\
* energy_overlap * scaling_factor
# if sum_density is true sum the orbital interactions
if sum_interaction == True:
new_overlap_orbitals = []
new_orbital_interaction = []
for count, orbital in enumerate(overlap_orbitals):
if orbital[0] not in new_overlap_orbitals:
new_overlap_orbitals.append(orbital[0])
new_orbital_interaction.append(orbital_interaction[count])
else:
new_orbital_interaction[-1] += orbital_interaction[count]
overlap_orbitals = np.array(new_orbital_interaction)
orbital_interaction = np.array(new_orbital_interaction)
return orbital_interaction
[docs] def _assign_orbitals(self, orbital_scores):
""" Assign gas orbitals to one or more adsorbate orbitals
Parameters
----------
orbital_scores : numpy.ndarray
M x N 2D array where M is the number of gas molecular orbitals
and N is the number of adsorbate molecular orbitals
Returns
-------
gas_2_adsorbate : numpy.ndarray
Array of gas and adsorbate orbital indices along with corresponding
band centers
"""
gas_band_centers = self.gas_band_centers
adsorbate_band_centers = self.adsorbate_band_centers
if orbital_scores.shape[0] >= orbital_scores.shape[1]:
orbital_assignments = np.argmax(orbital_scores,axis=1).astype('int')
indices = np.arange(orbital_assignments.size)
gas_2_adsorbate = np.vstack((indices,orbital_assignments\
, gas_band_centers, adsorbate_band_centers[orbital_assignments])).T
else:
orbital_assignments = np.argmax(orbital_scores,axis=0).astype('int')
indices = np.arange(orbital_assignments.size)
gas_2_adsorbate = np.vstack((orbital_assignments,indices\
, gas_band_centers[orbital_assignments], adsorbate_band_centers)).T
return gas_2_adsorbate
[docs] def get_orbital_scores(self, gas_features, adsorbate_features, energy_weight):
""" Select the set of orbital scores to be used
Parameters
----------
gas_orbital_features : numpy.ndarray
gas orbital features
adsorbate_orbital_features : numpy.ndarray
adsorbate orbital features
energy_weight : float
The degree to weight the energy error before computing the
associated liklihood that goes into computing the scores. Increasing
the energy weight sharpens the univariate distribution around
matching peak positions.
Returns
-------
orbital_scores : numpy.ndarray
M x N 2D array where M is the number of gas molecular orbitals
and N is the number of adsorbate molecular orbitals
"""
orbital_scores_gas\
= self._get_orbital_scores(gas_features, adsorbate_features, energy_weight)
orbital_scores_adsorbate\
= self._get_orbital_scores(adsorbate_features, gas_features, energy_weight)
if orbital_scores_gas.shape[0] == orbital_scores_adsorbate.shape[0]:
orbital_scores = (orbital_scores_gas * orbital_scores_adsorbate.T)**0.5
orbital_scores = (orbital_scores * orbital_scores.T)**0.5
elif orbital_scores_gas.shape[0] > orbital_scores_adsorbate.shape[0]:
orbital_scores = orbital_scores_gas
else:
orbital_scores = orbital_scores_adsorbate.T
return orbital_scores
[docs] @staticmethod
def _get_sum_score(orbital_scores):
""" Compute sum for optimizing upshift
Parameters
----------
orbital_scores : numpy.ndarray
M x N 2D array where M is the number of gas molecular orbitals
and N is the number of adsorbate molecular orbitals
Returns
-------
sum_score : numpy.ndarray
Sum of highest scores for each gas or adsorbate molecular orbital
"""
if orbital_scores.shape[0] >= orbital_scores.shape[1]:
sum_score = orbital_scores.max(axis=1).sum()
else:
sum_score = orbital_scores.max(axis=0).sum()
return sum_score
[docs] def _get_upshift_score(self, upshift):
""" Obtain orbital scores given an upshift value
Parameters
----------
upshift : float
The energy the gas molecular orbitals are shifted before pairing
Returns
-------
sum_score : numpy.ndarray
Sum of highest scores for each gas or adsorbate molecular orbital
"""
gas_features\
= self._get_orbital_features(self.GAS_PDOS\
, self.gas_orbital_indices\
, self.gas_indices, upshift)
adsorbate_features\
= self._get_orbital_features(self.REFERENCE_PDOS\
, self.adsorbate_orbital_indices\
, self.adsorbate_indices, 0)
scores = self.get_orbital_scores(gas_features, adsorbate_features\
, self.energy_weight)
sum_score = self._get_sum_score(scores)
return sum_score
[docs] def optimize_energy_shift(self, bound=(-0.5,1.5), reset=False, plot=False):
""" Get the optimal energy shift in the gas and adsorbate orbitals
Parameters
----------
upshift : float
The energy the gas molecular orbitals are shifted before pairing
bound : tuple
Lower and upper bounds to for which to find the optimal energy shift
rest : bool
Indicates whether the values contained in the PDOS_OVERLAP object
should be reset using the optimal energy shift
plot : bool
Indicates whether optimization results should be plotted
Returns
-------
optimized_upshift : float
The fraction of the fermi energy the gas and adsorbate molecular
orbitals are shifted before pairing for maximizing likelihood of
score pairings
"""
def f(upshift):
return -1 * self._get_upshift_score(upshift)
grid = np.linspace(bound[0],bound[1],50,endpoint=True)
y = np.array([f(x) for x in grid])
argmin = np.argmin(y)
opt = minimize_scalar(f, bounds=(grid[argmin - 2], grid[argmin + 2])\
, method='bounded', options={'xatol': 1e-05\
, 'maxiter': 500, 'disp': 1})
optimized_upshift = opt['x']
if reset == True:
self.__init__(self.GAS_PDOS, self.REFERENCE_PDOS\
, self.adsorbate_indices, self.site_indices\
, self.min_occupation, upshift=optimized_upshift\
, energy_weight=self.energy_weight\
, sum_density=False, sum_spin=True)
if plot == True:
set_figure_settings('paper')
plt.figure()
plt.plot(grid, - 1 * y)
plt.xlabel(r'Shift in gas molecular orbitals [eV]')
plt.ylabel(r'$\sum_{m=1}^{M}$max(orbital score$_{m,j}$ for j in N )')
plt.show()
return optimized_upshift
[docs] def plot_energy_overlap(self, indices = [...], atomic_orbitals=['s','p','d']
, sum_overlap=False, figure_directory='print',extension='jpg'):
""" Plot energy overlap of atomic and molecular orbitals
Parameters
----------
indices : list[int]
Index values of adsorbate molecular orbitals for which to plot
energy overlap histograms.
figure_directory : str
indicates how to display or save the figures or directory to save
atomic_orbitals : list[str]
Adsorption site atomic orbitals of interest.
extension : str
'pdf' or 'jpg' how to save the file
"""
get_energy_overlap = self.get_energy_overlap
overlap_orbitals, energy_overlap = get_energy_overlap(atomic_orbitals
, sum_overlap=sum_overlap
, sum_spin=self.sum_spin
, normalize=False)
energy_overlap = energy_overlap[indices]
if figure_directory not in ['presentation', 'paper', 'print']:
if len(indices) == 2:
fig = plt.figure(figsize=(7.2,2),dpi=400)
abc = ['(a)','(b)']
axes = fig.subplots(nrows=1, ncols=2)
axes_list = [axes[0], axes[1], axes[1,0]]
else:
fig = plt.figure(figsize=(7.2,4),dpi=400)
abc = ['(a)','(b)','(c)','(d)']
axes = fig.subplots(nrows=2, ncols=2)
axes_list = [axes[0,0], axes[0,1], axes[1,0], axes[1,1]]
xticks = np.arange(1,len(overlap_orbitals)+1)
#plotting function
for index, overlap in enumerate(energy_overlap):
axes_list[index].text(0.93,0.90,abc[index],transform=axes_list[index].transAxes)
axes_list[index].bar(xticks,overlap)
axes_list[index].set_xticks(xticks)
axes_list[index].set_xticklabels(overlap_orbitals)
fig.text(0.001, 0.5, 'Overlap with adsorbate molecular orbitals [states]', va='center', rotation='vertical')
fig.text(0.5, 0.01, 'Metal atomic orbitals', ha='center')
fig.set_tight_layout({'pad':2,'w_pad':1,'h_pad':1})
figure_path = os.path.join(figure_directory,'energy_overlap.'+extension)
plt.savefig(figure_path, format=extension)
plt.close()
else:
if figure_directory == 'presentation':
set_figure_settings('presentation')
if len(energy_overlap.shape) == 1:
energy_overlap = [energy_overlap]
for overlap in energy_overlap:
plt.figure()
xticks = np.arange(1,len(overlap_orbitals)+1)
plt.bar(xticks,overlap)
plt.xticks(xticks,overlap_orbitals)
plt.xlabel('Metal states projected onto atomic orbitals')
plt.ylabel('Overlap with adsorbate molecular orbitals')
plt.show()
[docs] def plot_projected_density(self,sum_density=True, sum_spin=True\
, figure_directory='print',extension='jpg'):
""" Plot projected density of the gas, adsorbate and site
Parameters
----------
sum_density : bool
Tag indicates if state density on the same sublevel should be
summed when calculating overlap of the site atomic orbitals and the
adsorbate molecular orbitals.
sum_spin : bool
Tag indictes if state density with different spins should be summed
when generating gas and adsorbate orbital features.
figure_directory : str
indicates how to display or save the figures
extension : str
'pdf' or 'jpg' how to save the file
"""
GAS_PDOS = self.GAS_PDOS
gas_indices = self.gas_indices
REFERENCE_PDOS = self.REFERENCE_PDOS
adsorbate_indices = self.adsorbate_indices
site_indices = self.site_indices
if figure_directory == 'presentation':
set_figure_settings('presentation')
else:
set_figure_settings('paper')
orbital_list=['s', 'p', 'd']
colors = ['b-','g-','r-']
zorder = [2,3,4]
if sum_spin == False:
orbital_list=['s+', 's-', 'p+', 'p-', 'd+', 'd-']
colors = ['b-', 'b-', 'g-', 'g-', 'r-', 'r-']
zorder = [2,3,4,5,6,7]
if figure_directory not in ['print', 'presentation']:
fig = plt.figure(figsize=(7.2,4),dpi=400)
else:
fig = plt.figure()
abc = ['(a)','(b)','(c)']
axes = fig.subplots(nrows=1, ncols=3)
axes_list = [axes[0], axes[1], axes[2]]
#plotting function
def plot_density(PDOS, indices, index):
orbitals, projected_density = PDOS.get_site_dos(indices\
, orbital_list, sum_density, sum_spin)
if sum_spin == False:
projected_density[1::2,:] *= -1
for count, density in enumerate(projected_density):
axes_list[index].plot(density, PDOS.get_energies(), colors[count], zorder=zorder[count])
axes_list[index].plot([np.min(projected_density), np.max(projected_density)]\
,[PDOS.e_fermi, PDOS.e_fermi]\
,'k--', zorder=1,linewidth=5)
#axes_list[index].legend([i for i in orbitals]+ ['fermi level'])
axes_list[index].text(0.90,0.96,abc[index],transform=axes_list[index].transAxes)
#plot gas density
plot_density(GAS_PDOS, gas_indices, 0)
#plot adsorbate density
plot_density(REFERENCE_PDOS, adsorbate_indices, 1)
#plot adsorption-site density
plot_density(REFERENCE_PDOS, site_indices, 2)
fig.text(0.001, 0.5, 'Energy [eV]', va='center', rotation='vertical')
fig.text(0.5, 0.01, 'State density [states/eV]', ha='center')
if figure_directory not in ['print', 'presentation']:
figure_path = os.path.join(figure_directory,'pdos.'+extension)
fig.set_tight_layout({'pad':2,'w_pad':1})
plt.savefig(figure_path, format=extension)
plt.close()
else:
plt.show()