# -*- coding: utf-8 -*-
import numpy
from params import tb_parm
from params import consts
from lib import libtb
TRUE_NAME = {"Carbon" : "C", "Hydrogen" : "H"}
class TBException(Exception):
def __init__(self, message):
super(TBException, self).__init__(message)
[docs]class TB(object):
"""
Class that realizes Tight-Binding method.
Constructor of class TB.
Args:
struct (Structure): nanostructure.
lib (module): dynamic library, it realize filling Hamiltonian matrix, calculating repulsive energie, forces acting on atoms.
parms (dict): parameters of method.
Kwargs:
functions (list): list of functions that will be loaded from lib library.
"""
def __init__(self, struct, parm):
if parm:
self.parm = tb_parm.getTB_parm(parm[0])
else:
self.parm = tb_parm.getTB_parm("YW")
functions = ["fill_matrix", "e_rep"]
self.struct = struct
if self.preprocessing(struct):
self.load_tb_functions(functions, libtb)
self.load_force_function(libtb)
else:
raise TBException("Preprocessing failed.")
[docs] def preprocessing(self, struct):
"""
Preprocessing for TB method, that:
#. Checks existing of parameters for different atom types, involved in structure.
#. Creates Hamiltonian matrix.
#. Calculates number of orbitals of each type and shift indices for orbitals in hamiltonian.
#. Calculates number of occupied energy levels
Args:
struct (Structure): nanostructure
"""
atypes = self.parm["atypes"]
if not hasattr(struct, "atom_types"):
struct.atom_types = numpy.zeros(struct.n_atoms, numpy.int32)
calc_mass = False
if not hasattr(struct, "mass"):
calc_mass = True
struct.mass = numpy.zeros(struct.n_atoms, numpy.float)
for i, name in enumerate(struct.atom_names):
if name in atypes:
struct.atom_types[i] = atypes.index(name)
else:
if name in TRUE_NAME:
struct.atom_types[i] = atypes.index(TRUE_NAME[name])
if calc_mass:
struct.mass[i] = consts.ATOM_MASS[TRUE_NAME[name]]
else:
raise TBException("Wrong atom type name %s" % name)
n_orbs = self.parm["n_orbs"]
struct.orb_cnt = n_orbs[struct.atom_types]
rank = struct.orb_cnt.sum()
struct.hami = numpy.zeros((rank, rank), numpy.float)
struct.rank = int(rank)
# TODO real value of n_occup
struct.n_occup = int(rank / 2 + rank % 2)
struct.n_atoms = int(struct.n_atoms)
struct.orb_inds = numpy.zeros(struct.orb_cnt.max(), numpy.int32)
struct.tb_links = struct.full_links
struct.num_tb_links = struct.num_full_links
for orbi in struct.orb_cnt:
struct.orb_inds[:orbi] += 1
struct.orb_inds[1:] += struct.orb_inds[:-1]
struct.orb_inds[1:] = struct.orb_inds[:-1]
struct.orb_inds[0] = 0
return True
[docs] def update_tb_links(self, tb_links, num_tb_links):
"""
#. Updates data of tight bonds in structure
#. Zeros hamiltonian matrix
Args:
tb_links(numpy.array): 2D arrays of bonds.
num_tb_links(numpy.array): 1D array, it contains numbers of bonds for each atom
"""
self.struct.tb_links = tb_links
self.struct.num_tb_links
self.hami.fill(0)
[docs] def load_tb_functions(self, functions, lib):
"""
Loading of functions from dynamic library and creating wrap-methods.
Args:
functions (list): list of functions, that will be loaded from lib
lib (module): dynamic library
"""
def get_func(name, lib):
def func():
return getattr(lib, name)(self.struct.__dict__, self.parm)
return func
for f_type in functions:
setattr(self, f_type, get_func(f_type, lib))
[docs] def load_force_function(self, lib):
"""
Loading of function that calculate forces and creating wrap-methods from dynamic library.
Args:
functions (list): list of functions, that will be loaded from lib
lib (module): dynamic library
"""
def get_func(name, lib):
def func(f, a, b):
return getattr(lib, name)(self.struct.__dict__, self.parm, f, a, b)
return func
if hasattr(lib, "forces"):
func = get_func("forces", lib)
setattr(self, "forces", func)
def update_atom_charges(self):
if not hasattr(self.struct, "eigvec"):
self.compute_ebond()
q = self.calc_orbital_charges(self.struct.eigvec)
self.struct.charges = numpy.zeros(self.struct.n_atoms, numpy.float)
for atom_id, at in enumerate(self.struct.atom_types):
atom_orb_count = self.struct.orb_cnt[atom_id]
atom_orb_indexes = [atom_id + self.struct.orb_inds[orb] for orb in range(atom_orb_count)]
self.struct.charges[atom_id] = q[atom_orb_indexes].sum() - atom_orb_count
def calc_orbital_charges(self, eig_vects):
return numpy.sum(eig_vects[:,:self.struct.n_occup] ** 2, axis = 1) * 2
[docs] def compute_erep(self):
r"""
Repulsive energy between two atoms
.. math::
:nowrap:
$E_{core}(r) = E_{core}(r_0) \left(\frac{r_0}{r}\right)^{m_a}exp\left(-m_b\left(\frac{r}{r_c}\right)^{m_c} + m_b\left( \frac{r_0}{r_c}\right)^{m_c}\right)$
Total repulsive energy
.. math::
:nowrap:
$E_{rep} = \sum\limits_{i < j}{E_{core}(r_{ij})}$
Returns:
energy_rep(float): repulsive energy
"""
return self.e_rep()
[docs] def compute_ebond(self):
r"""
Calculating bond energy using formula
.. math::
:nowrap:
$E_{val} = \sum\limits_{\gamma}{n_\gamma \epsilon_\gamma}$
where
:math:`\epsilon_{\gamma}`
is a value on diagonal of diagonalised hamiltonian
Returns:
energy_bond(float): bond energy of atoms
"""
self.fill_matrix()
self.struct.eigval, self.struct.eigvec = numpy.linalg.eigh(self.struct.hami)
e = 2. * self.struct.eigval[:self.struct.n_occup].sum()
return e
[docs] def compute(self):
"""
Calculating full energy of nanostructure.
Returns:
energy_total(float): full energy of nanostructure
"""
return self.compute_ebond() + self.compute_erep()
[docs] def ensure_forces(self):
"""
Function changes method for calculating forces, acting on atoms of
structure, to ensure, that hamiltonian matrix is not zero before
calculation start and eigenvectors can be calculated.
"""
forces = self.forces
def wrapped(f, a, b):
self.compute_ebond()
forces(f, a, b)
self.forces = wrapped