Source code for kvazar.core.tb

# -*- 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 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