From 287663dbb71d9c9712e1bbc1975e45e2ccfffbde Mon Sep 17 00:00:00 2001 From: Jijie Zou <73353910+AsymmetryChou@users.noreply.github.com> Date: Mon, 11 Dec 2023 16:17:38 +0800 Subject: [PATCH] add TBtrans Interface tbtrans_init.py and its test files (#50) * add tbtrans_init and related files for Shell command * rename load_dptb_model as load_model * add or modify some docstring * add test files for tbtrans_init.py * modify test_tbtrans_init.py * add some comments * change test_hBN_struct to orthogonal unit cell * remove unnecessary sisl import * add warning for input structure file in vasp format. * add docstring to class * fix some small problems * rewrite orbitals_get function * hamil_get to hamil_get_write * add runtime error for .vasp and change R_vec assert * add comments to sort([2,1,0]) * add TODO for double leads and notes on Z-direction trans. * move set_nsc to hamiltonian_get * add notes on Hamil_sisl.set_nsc * add import properties in __init__ * remove unnecessary print * rewrite load_model in for * fix load_model docstring problem * delete load_model and merge its function to hamil_get_write * delete load_model * fix bugs in chemical_symbol determination * check _orbital_name_get --- dptb/entrypoints/run.py | 9 + dptb/postprocess/tbtrans_init.py | 664 ++++++++++++++++++ dptb/tests/data/test_negf/lead_R.nc | Bin 0 -> 61494 bytes .../best_nnsk_b3.600_c3.600_w0.300.json | 122 ++++ .../tests/data/test_tbtrans/negf_tbtrans.json | 90 +++ .../test_tbtrans/test_hBN_zigzag_struct.xyz | 22 + dptb/tests/test_tbtrans_init.py | 106 +++ dptb/utils/argcheck.py | 58 +- 8 files changed, 1070 insertions(+), 1 deletion(-) create mode 100644 dptb/postprocess/tbtrans_init.py create mode 100644 dptb/tests/data/test_negf/lead_R.nc create mode 100644 dptb/tests/data/test_tbtrans/best_nnsk_b3.600_c3.600_w0.300.json create mode 100644 dptb/tests/data/test_tbtrans/negf_tbtrans.json create mode 100644 dptb/tests/data/test_tbtrans/test_hBN_zigzag_struct.xyz create mode 100644 dptb/tests/test_tbtrans_init.py diff --git a/dptb/entrypoints/run.py b/dptb/entrypoints/run.py index c9ba3a83..afc7ccf5 100644 --- a/dptb/entrypoints/run.py +++ b/dptb/entrypoints/run.py @@ -23,6 +23,7 @@ from dptb.postprocess.bandstructure.ifermi_api import ifermiapi, ifermi_installed, pymatgen_installed from dptb.postprocess.write_skparam import WriteNNSKParam from dptb.postprocess.NEGF import NEGF +from dptb.postprocess.tbtrans_init import TBTransInputSet,sisl_installed __all__ = ["run"] @@ -215,6 +216,14 @@ def run( negf.compute() log.info(msg='NEGF calculation successfully completed.') + if task == 'tbtrans_negf': + if not(sisl_installed): + log.error(msg="sisl is required to perform tbtrans calculation !") + raise RuntimeError + + tbtrans_init = TBTransInputSet(apiHrk, run_opt, task_options) + tbtrans_init.hamil_get_write(write_nc=True) + log.info(msg='TBtrans input files are successfully generated.') if output: with open(os.path.join(output, "run_config.json"), "w") as fp: diff --git a/dptb/postprocess/tbtrans_init.py b/dptb/postprocess/tbtrans_init.py new file mode 100644 index 00000000..8814b071 --- /dev/null +++ b/dptb/postprocess/tbtrans_init.py @@ -0,0 +1,664 @@ +import numpy as np +import logging +import shutil +import re +from dptb.structure.structure import BaseStruct +from dptb.utils.tools import j_loader,j_must_have + +from ase.io import read,write +from ase.build import sort +import ase.atoms +import torch + +from dptb.utils.constants import atomic_num_dict_r + + +log = logging.getLogger(__name__) + +try: + import sisl + from sisl.orbital import Orbital + + sisl_installed = True + +except ImportError: + log.error('sisl is not installed.Thus the input for TBtrans can not be generated, please install it first!') + sisl_installed = False + + +if shutil.which('tbtrans') is None: + log.error('tbtrans is not in the Environment PATH. Thus the input for TBtrans can be generated but not run.') + +# TBTransInputSet is used to transform input data for DeePTB-negf into TBtrans input files. +# TBtrans (Tight-Binding transport) is a generic computer program which calculates transport and other physical quantities +# using the Green function formalism. It is a stand-alone program which allows extreme scale tight-binding calculations. +# For details, see https://www.sciencedirect.com/science/article/pii/S001046551630306X?via%3Dihub. +# To run TBTransInputSet, user need sisl package(https://zerothi.github.io/sisl/index.html) + + +class TBTransInputSet(object): + """ The TBTransInputSet class is used to transform input data for DeePTB-negf into a TBTrans object. + + Attention: the transport direction is forced to be z direction in this stage, please make sure the structure is in + correct direction. + + Properties + ----------- + - apiHrk + apiHrk has been loaded in the run.py file. It is used as an API for + performing certain operations or accessing certain functionalities. + - run_opt + The `run_opt` parameter is a dictionary that contains options for running the model. + It has been loaded and prepared in the run.py file. + - jdata + jdata is a JSON object that contains options and parameters for the task Generation of Input Files for TBtrans. + It is loaded in the run.py. + - results_path + The `results_path` parameter is a string that represents the path to the directory where the + results will be saved. + - stru_options + The `stru_options` parameter is a dictionary that contains options for the structure from DeePTB input. + - energy_unit_option + The `energy_unit_option` parameter is a string that specifies the unit of energy for the + calculation. It can be either "Hartree" or "eV". + - geom_all + The `geom_all` parameter is the geometry of the whole structure, including the device and leads. + - H_all + The `H_all` parameter is the sisl.Hamiltonian for the entire system, including the device and leads. + - H_lead_L + The `H_lead_L` parameter is sisl.Hamiltonian for the left lead. + - H_lead_R + The `H_lead_R` parameter is sisl.Hamiltonian for the right lead. + - allbonds_all + The `allbonds_all` parameter is a tensor that contains all of the bond information for the entire system. + - allbonds_lead_L + The `allbonds_lead_L` parameter is a tensor that contains all of the bond information for the left lead. + - allbonds_lead_R + The `allbonds_lead_R` parameter is a tensor that contains all of the bond information for the right lead. + - hamil_block_all + The `hamil_block_all` parameter is a tensor that contains the Hamiltonian matrix elements for each specific bond in allbonds_all. + - hamil_block_lead_L + The `hamil_block_lead_L` parameter is a tensor that contains the Hamiltonian matrix elements for each specific bond in allbonds_lead_L. + - hamil_block_lead_R + The `hamil_block_lead_L` parameter is a tensor that contains the Hamiltonian matrix elements for each specific bond in allbonds_lead_R. + - overlap_block_all + The `overlap_block_all` parameter is a tensor that contains the overlap matrix elements for each specific basis in the entire system. + - overlap_block_lead_L + The `overlap_block_lead_L` parameter is a tensor that contains the overlap matrix elements for each specific basis in the left lead. + - overlap_block_lead_R + The `overlap_block_lead_R` parameter is a tensor that contains the overlap matrix elements for each specific basis in the right lead. + """ + def __init__(self, apiHrk, run_opt, jdata): + self.apiHrk = apiHrk #apiHrk has been loaded in run.py + self.jdata = jdata #jdata has been loaded in run.py, jdata is written in negf.json + + self.results_path = run_opt['results_path'] + if not self.results_path.endswith('/'):self.results_path += '/' + self.stru_options = j_must_have(jdata, "stru_options") + self.energy_unit_option = 'eV' # enenrgy unit for TBtrans calculation + + + self.geom_all,self.geom_lead_L,self.geom_lead_R,self.all_tbtrans_stru,self.lead_L_tbtrans_stru,self.lead_R_tbtrans_stru\ + = self.read_rewrite_structure(run_opt['structure'],self.stru_options,self.results_path) + + + self.orbitals_get(self.geom_all,self.geom_lead_L,self.geom_lead_R,apiHrk=apiHrk) + + self.H_all = sisl.Hamiltonian(self.geom_all) + self.H_lead_L = sisl.Hamiltonian(self.geom_lead_L) + self.H_lead_R = sisl.Hamiltonian(self.geom_lead_R) + + + #important properties for later use + + ##allbonds matrx, hamiltonian matrix, overlap matrix for the whole structure + self.allbonds_all = None + self.hamil_block_all = None + self.overlap_block_all = None + ##allbonds matrx, hamiltonian matrix, overlap matrix for lead_L + self.allbonds_lead_L = None + self.hamil_block_lead_L = None + self.overlap_block_lead_L = None + ##allbonds matrx, hamiltonian matrix, overlap matrix for lead_R + self.allbonds_lead_R = None + self.hamil_block_lead_R = None + self.overlap_block_lead_R = None + + + def hamil_get_write(self,write_nc:bool=True): + + '''The function `hamil_get_write` loads models for different structure.retrieves the Hamiltonian and overlap matrices / + for the device and left and right leads, then writes the contents of `self.H_all`, `self.H_lead_L`, and `self.H_lead_R` to nc files for + TBtrans calculations. + + `all` refers to the entire system, including the device and leads. + + + Returns + ------- + - allbonds_all: all of the bond information + - hamil_block_all: Hamiltonian block for the entire system, which is a tensor that contains + the values of the Hamiltonian matrix elements for each specific bond in allbonds_all + - overlap_block_all: overlap block for the entire system, which is a tensor that contains + the values of the overlap matrix elements for each specific basis + ''' + + + # get the Hamiltonian matrix for the entire system + self.allbonds_all,self.hamil_block_all,self.overlap_block_all\ + =self._load_model(self.apiHrk,self.all_tbtrans_stru) + self.hamiltonian_get(self.allbonds_all,self.hamil_block_all,self.overlap_block_all,\ + self.H_all,self.energy_unit_option) + + # get the Hamiltonian matrix for the left lead + self.allbonds_lead_L,self.hamil_block_lead_L,self.overlap_block_lead_L\ + =self._load_model(self.apiHrk,self.lead_L_tbtrans_stru) + self.hamiltonian_get(self.allbonds_lead_L,self.hamil_block_lead_L,self.overlap_block_lead_L,\ + self.H_lead_L,self.energy_unit_option) + + # get the Hamiltonian matrix for the right lead + self.allbonds_lead_R,self.hamil_block_lead_R,self.overlap_block_lead_R\ + =self._load_model(self.apiHrk,self.lead_R_tbtrans_stru) + self.hamiltonian_get(self.allbonds_lead_R,self.hamil_block_lead_R,self.overlap_block_lead_R,\ + self.H_lead_R,self.energy_unit_option) + + if write_nc: + self.H_all.write(self.results_path+'structure.nc') + self.H_lead_L.write(self.results_path+'lead_L.nc') + self.H_lead_L.write(self.results_path+'lead_R.nc') + else: + print('Hamiltonian matrices have been generated, but not written to nc files(TBtrans input file).') + + # def hamil_write(self): + # '''The function `hamil_write` writes the contents of `self.H_all`, `self.H_lead_L`, and `self.H_lead_R` + # to separate files in the `results_path` directory. + + # ''' + # self.H_all.write(self.results_path+'structure.nc') + # self.H_lead_L.write(self.results_path+'lead_L.nc') + # self.H_lead_L.write(self.results_path+'lead_R.nc') + + + + def read_rewrite_structure(self,structure_file:str,struct_options:dict,results_path:str): + '''The function `read_rewrite_structure` reads a structure file, extracts specific regions of the structure, + sorts the atoms in the structure, and outputs the sorted structures in XYZ and VASP file formats for later operations. + + Parameters + ---------- + structure_file + The `structure_file` parameter is the path to the file containing the structure information of the + system you want to analyze. It can be in either VASP format (.vasp) or XYZ format (.xyz). + struct_options + The `struct_options` parameter is a dictionary that contains various options for the structure. + result_path + The `result_path` parameter is the path where the output files will be saved. + + Returns + ------- + + - geom_all: the geometry of the entire structure,including the device and leads + - geom_lead_L: the geometry of the left lead + - geom_lead_R: the geometry of the right lead + - all_tbtrans_stru: the path to the sorted xyz file for the entire structure + - lead_L_tbtrans_stru: the path to the sorted xyz file for the left lead + - lead_R_tbtrans_stru: the path to the sorted xyz file for the right lead + + ''' + + + lead_L_id=struct_options['lead_L']['id'].split('-') + lead_R_id=struct_options['lead_R']['id'].split('-') + # device_id=struct_options['device']['id'].split('-') + + lead_L_range=[i for i in range(int(lead_L_id[0]), int(lead_L_id[1]))] + lead_R_range=[i for i in range(int(lead_R_id[0]), int(lead_R_id[1]))] + + + # Structure input: read vasp file + # structure_vasp = sisl.io.carSileVASP(structure_file) + # geom_device = structure_vasp.read_geometry() + if structure_file.split('.')[-1]=='vasp': + structure_vasp = sisl.io.carSileVASP(structure_file) + geom_all = structure_vasp.read_geometry() + if geom_all.atoms.nspecie>1: + raise RuntimeError('ERROR! In transport calculation, VASP structure file is only valid for materials with one single element!') + elif structure_file.split('.')[-1]=='xyz': + structure_xyz = sisl.io.xyzSile(structure_file) + geom_all = structure_xyz.read_geometry() + else: + raise RuntimeError('Structure file format is not supported. Only support vasp and xyz format') + # structure_xyz = sisl.io.xyzSile(structure_file) + # geom_device = structure_xyz.read_geometry() + #define lead geometry structure + geom_lead_L = geom_all.sub(lead_L_range) # subset of the geometry + geom_lead_R = geom_all.sub(lead_R_range) + + #sort sturcture atoms according to z-direction: it's easier for later coding + #2,1,0 refers to sort by axis=0 firstly, then axis=1, last for axis=2 + geom_lead_R = geom_lead_R.sort(axis=(2,1,0));geom_lead_L=geom_lead_L.sort(axis=(2,1,0)) + geom_all=geom_all.sort(axis=(2,1,0)) + + ##redefine the Lattice vector of Lead L/R + # lead_L_cor = geom_lead_L.axyz() + # Natom_PL = int(len(lead_L_cor)/2) + # first_PL_leadL = lead_L_cor[Natom_PL:];second_PL_leadL =lead_L_cor[:Natom_PL] + # PL_leadL_zspace = first_PL_leadL[0][2]-second_PL_leadL[-1][2] # the distance between Principal layers + # geom_lead_L.lattice.cell[2,2]=first_PL_leadL[-1][2]-second_PL_leadL[0][2]+PL_leadL_zspace + # assert geom_lead_L.lattice.cell[2,2]>0 + + lead_L_cor = geom_lead_L.axyz() #Return the atomic coordinates in the supercell of a given atom. + cell = np.array(geom_lead_L.lattice.cell)[:2] + Natom_PL = int(len(lead_L_cor)/2) + first_PL_leadL = lead_L_cor[Natom_PL:];second_PL_leadL =lead_L_cor[:Natom_PL] + R_vec = first_PL_leadL - second_PL_leadL + # assert np.abs(R_vec[0] - R_vec[-1]).sum() < 1e-5 + assert np.abs(R_vec[0] - R_vec.mean(axis=0)).sum() < 1e-5 + R_vec = R_vec.mean(axis=0) * 2 + cell = np.concatenate([cell, R_vec.reshape(1,-1)]) + # PL_leadL_zspace = first_PL_leadL[0][2]-second_PL_leadL[-1][2] # the distance between Principal layers + geom_lead_L.lattice.cell=cell + # assert geom_lead_L.lattice.cell[2,2]>0 + + #TODO: This version tbtrans_init only supports double lead case. Code for more leads will be added later. + + lead_R_cor = geom_lead_R.axyz() + cell = np.array(geom_lead_R.lattice.cell)[:2] + Natom_PL = int(len(lead_R_cor)/2) + first_PL_leadR = lead_R_cor[:Natom_PL];second_PL_leadR = lead_R_cor[Natom_PL:] + R_vec = first_PL_leadR - second_PL_leadR + assert np.abs(R_vec[0] - R_vec.mean(axis=0)).sum() < 1e-5 + R_vec = -1*R_vec.mean(axis=0) * 2 + cell = np.concatenate([cell, R_vec.reshape(1,-1)]) + # PL_leadR_zspace = second_PL_leadR[0][2]-first_PL_leadR[-1][2] + geom_lead_R.lattice.cell = cell + # print(cell) + # assert geom_lead_R.lattice.cell[2,2]>0 + + # set supercell + # PBC requirement in TBtrans + ## lead calculation have periodicity in all directions,which is different from dptb-negf + ## all(lead + central part) have periodicity in x,y,z directions: interaction between supercells + ### not sure that geom_all need pbc in z direction + + # pbc = struct_options['pbc'] + # if pbc[0]==True: nsc_x = 3 + # else: nsc_x = 1 + + # if pbc[1]==True: nsc_y = 3 + # else: nsc_y = 1 + + # geom_lead_L.set_nsc(a=nsc_x,b=nsc_y,c=3) #Set the number of super-cells in the `Lattice` object + # geom_lead_R.set_nsc(a=nsc_x,b=nsc_y,c=3) + # geom_all.set_nsc(a=nsc_x,b=nsc_y,c=3) + + + # output sorted geometry into xyz Structure file + all_tbtrans_stru=results_path+'structure_tbtrans.xyz' + sorted_structure = sisl.io.xyzSile(all_tbtrans_stru,'w') + geom_all.write(sorted_structure) + + lead_L_tbtrans_stru=results_path+'lead_L_tbtrans.xyz' + sorted_lead_L = sisl.io.xyzSile(lead_L_tbtrans_stru,'w') + geom_lead_L.write(sorted_lead_L) + + lead_R_tbtrans_stru=results_path+'lead_R_tbtrans.xyz' + sorted_lead_R = sisl.io.xyzSile(lead_R_tbtrans_stru,'w') + geom_lead_R.write(sorted_lead_R) + + # output sorted geometry into vasp Structure file: rewrite xyz files + ## writen as VASP for VESTA view + all_struct = read(all_tbtrans_stru) + all_vasp_struct = results_path+'structure_tbtrans.vasp' + write(all_vasp_struct,sort(all_struct),format='vasp') + + return geom_all,geom_lead_L,geom_lead_R,all_tbtrans_stru,lead_L_tbtrans_stru,lead_R_tbtrans_stru + + + + def orbitals_get(self,geom_all, geom_lead_L,geom_lead_R,apiHrk): + '''The function `orbitals_get` takes in various inputs such as geometric devices, leads, deeptb model, and + configurations, and assigns orbitals number, orbital names, shell-electron numbers to the atoms in the given sisl geometries . + + Here the geometry class is sisl.geometry, which is different from the structure class in dptb-negf. + We initialize sisl.geometry from structure files directly, therefore there is no orbital information in sisl.geometry. + + Parameters + ---------- + geom_all + The `geom_all` parameter is the geometry of the whole structure. It contains information about the + atoms and their positions. + geom_lead_L + The `geom_lead_L` parameter represents the geometry of the left lead. It contains + information about the atoms and their positions in the lead. + geom_lead_R + The `geom_lead_R` parameter represents the geometry of the right lead. It contains + information about the atoms in the lead, such as their positions and chemical symbols. + apiHrk + apiHrk has been loaded in the run.py file. It is used as an API for + performing certain operations or accessing certain functionalities when loading dptb model. + + ''' + n_species_lead_L = geom_lead_L.atoms.nspecie + n_species_lead_R = geom_lead_R.atoms.nspecie + n_species_all = geom_all.atoms.nspecie + n_species_list = [n_species_lead_L,n_species_lead_R,n_species_all] + geom_list = [geom_lead_L,geom_lead_R,geom_all] + + + dict_element_orbital = apiHrk.apihost.model_config['proj_atom_anglr_m'] + dict_shell_electron = apiHrk.apihost.model_config['proj_atom_neles'] + + for n_species, geom_part in zip(n_species_list,geom_list): + # species_symbols = split_string(geom_part.atoms.formula()) + ## get the chemical symbol of the part + # species_symbols = ''.join(char for char in geom_part.atoms.formula() if char.isalpha()) + + uni_symbol_index = np.unique(geom_part.atoms.Z) + species_symbols=[atomic_num_dict_r[i] for i in uni_symbol_index] + + assert len(species_symbols)==n_species # number of chemical elements in this part + + for i in range(n_species): #determine the orbitals number for each species + element_orbital_list = dict_element_orbital[species_symbols[i]] + # Examples of elemet_orbital_list: ['3s', '3p', 'd*'] + element_orbital_name = self._orbitals_name_get(element_orbital_list) + # Examples of element_orbital_name: ['3s', '3py', '3pz', '3px', 'dxy*', 'dyz*', 'dz2*', 'dxz*', 'dx2-y2*'] + + # shell_elec_num = self._shell_electrons(species_symbols[i]) + shell_elec_num = dict_shell_electron[species_symbols[i]] + shell_elec_list = np.zeros(len(element_orbital_name)) + shell_elec_list[:shell_elec_num]=1 #occupation number for each orbital, here we assume the system is in ground state + + for atom_index in range(geom_part.na): + if geom_part.atoms[atom_index].symbol == species_symbols[i]: + geom_part.atoms[atom_index]._orbitals =[Orbital(-1, q0=q,tag=tag) for q,tag in zip(shell_elec_list,element_orbital_name)] + # attribute sisl.Orbital object to each atom in sisl.geometry.atoms[x]._orbitals + + geom_lead_L.atoms._update_orbitals() #sisl use ._update_orbitals() to ensure the order of orbitals + geom_lead_R.atoms._update_orbitals() + geom_all.atoms._update_orbitals() + + + def _orbitals_name_get(self,element_orbital_class:str): + '''The `_orbitals_name_get` function takes a list of element orbital classes and returns a list of + orbital names. + + Parameters + ---------- + element_orbital_class + A list of strings representing the orbital classes of an element. Each string in the list + represents a different orbital class. + + Returns + ------- + a list of orbital names. + + Examples + -------- + >>>element_orbital_name = self._orbitals_name_get(['2s', '2p']) + >>>element_orbital_name + ['2s', '2py', '2pz', '2px'] + >>>element_orbital_name = self._orbitals_name_get(['3s', '3p', 'd*']) + >>>element_orbital_name + ['3s', '3py', '3pz', '3px', 'dxy*', 'dyz*', 'dz2*', 'dxz*', 'dx2-y2*'] + ''' + orbital_name_list=[] + for orb_cla in element_orbital_class: # ['3s', '3p', 'd*'] + orbital_name = re.findall(r'\d+|[a-zA-Z]+|\*',orb_cla) + orbital_name = sorted(orbital_name, key=lambda x: (x == '*',x.isnumeric(),x))# put s,p,d in the 1st position + + assert len(orbital_name)>1 #example 3d*: orbital_name: ['d', '3', '*'] + if orbital_name[0]=='s': + orbital_name_list += [orbital_name[1]+'s'] + elif orbital_name[0]=='p': + if orbital_name[1].isnumeric(): # not polarized orbital + orbital_name_list += [orbital_name[1]+'py',orbital_name[1]+'pz',orbital_name[1]+'px'] + else:#polarized orbital + orbital_name_list += ['py*','pz*','px*'] + elif orbital_name[0]=='d': + if orbital_name[1].isnumeric(): # not polarized orbital + orbital_name_list += [orbital_name[1]+'dxy',orbital_name[1]+'dyz',\ + orbital_name[1]+'dz2',orbital_name[1]+'dxz',orbital_name[1]+'dx2-y2'] + else:#polarized orbital + orbital_name_list += ['dxy*','dyz*','dz2*','dxz*','dx2-y2*'] + else: + raise RuntimeError("At this stage dptb-negf only supports s, p, d orbitals") + # print(orbital_name_list) + # raise RuntimeError('stop here') + return orbital_name_list + + # def _shell_electrons(self,element_symbol): + # '''The function `_shell_electrons` calculates the number of shell electrons for a given element symbol. + + # In this code, shell electron number is trivial for subgroup element. It would be improved soon. + + # Parameters + # ---------- + # element_symbol + # The element symbol is a string representing the symbol of an element on the periodic table. For + # example, "H" for hydrogen, "O" for oxygen, or "Fe" for iron. + + # Returns + # ------- + # the number of shell electrons for the given element symbol. + + # ''' + # atomic_number = PeriodicTable().Z_int(element_symbol) + # assert atomic_number > 1 and atomic_number <=118 + + # if atomic_number>18: + # print('In this code, shell electron number is trivial for subgroup element ') + # rare_element_index = [2,10,18,36,54,86] + + # for index in range(len(rare_element_index)-1): + # if atomic_number > rare_element_index[index] and atomic_number <= rare_element_index[index+1]: + # core_ele_num = atomic_number-rare_element_index[index] + + # print(element_symbol+' shell elec: '+str(core_ele_num)) + # return core_ele_num + + + + # def _load_dptb_model(self,checkfile:str,config:str,structure_tbtrans_file:str,run_sk:bool,use_correction:Optional[str]): + + # def _load_model(self,apiHrk,structure_tbtrans_file:str): + # '''The `_load_model` function loads model from deeptb and returns the Hamiltonian elements. + + # Parameters + # ---------- + # apiHrk + # apiHrk has been loaded in the run.py file. It is used as an API for + # performing certain operations or accessing certain functionalities when loading dptb model. + # structure_tbtrans_file : str + # The parameter `structure_tbtrans_file` is a string that represents the file path to the structure + # file in the TBTrans format. + + # Returns + # ------- + # The function `_load_model` returns three variables: `allbonds`, `hamil_block`, and + # `overlap_block`. + + # ''' + # if all((use_correction, run_sk)): + # raise RuntimeError("--use-correction and --train_sk should not be set at the same time") + + # ## read Hamiltonian elements + # if run_sk: + # apihost = NNSKHost(checkpoint=checkfile, config=config) + # apihost.register_plugin(InitSKModel()) + # apihost.build() + # ## define nnHrk for Hamiltonian model. + # apiHrk = NN2HRK(apihost=apihost, mode='nnsk') + # else: + # apihost = DPTBHost(dptbmodel=checkfile,use_correction=use_correction) + # apihost.register_plugin(InitDPTBModel()) + # apihost.build() + # apiHrk = NN2HRK(apihost=apihost, mode='dptb') + + + # self.allbonds_all,self.hamil_block_all,self.overlap_block_all\ + # =self._load_model(self.apiHrk,self.all_tbtrans_stru) + # self.allbonds_lead_L,self.hamil_block_lead_L,self.overlap_block_lead_L\ + # =self._load_model(self.apiHrk,self.lead_L_tbtrans_stru) + # self.allbonds_lead_R,self.hamil_block_lead_R,self.overlap_block_lead_R\ + # =self._load_model(self.apiHrk,self.lead_R_tbtrans_stru) + + # structure_tbtrans_file_list = [self.all_tbtrans_stru,self.lead_L_tbtrans_stru,self.lead_R_tbtrans_stru] + + ## create BaseStruct + # structure_base =BaseStruct( + # atom=ase.io.read(structure_tbtrans_file), + # format='ase', + # cutoff=apiHrk.apihost.model_config['bond_cutoff'], + # proj_atom_anglr_m=apiHrk.apihost.model_config['proj_atom_anglr_m'], + # proj_atom_neles=apiHrk.apihost.model_config['proj_atom_neles'], + # onsitemode=apiHrk.apihost.model_config['onsitemode'], + # time_symm=apiHrk.apihost.model_config['time_symm'] + # ) + + # apiHrk.update_struct(structure_base) + # allbonds,hamil_block,overlap_block = apiHrk.get_HR() + + # return allbonds,hamil_block,overlap_block + + + + + def _load_model(self,apiHrk,structure_tbtrans_file:str): + '''The `load_dptb_model` function loads a DPTB or NNSK model and returns the Hamiltonian elements. + Here run_sk is a boolean flag that determines whether to run the model using the NNSK or DPTB. + + Parameters + ---------- + checkfile : str + The `checkfile` parameter is the file path to the model checkpoint file. + config : str + The `config` parameter is a string that represents the configuration file for the model. It + contains information such as the model architecture, hyperparameters, and other settings that are + necessary for loading and building the model. + structure_tbtrans_file : str + The `structure_tbtrans_file` parameter is the file path to the structure file in the TBtrans + format. + struct_option : dict + The `struct_option` parameter is a dictionary that contains various options for the structure. It + includes the following keys: + run_sk : bool + The `run_sk` parameter is a boolean flag that determines whether to run the model using the NNSK + (Neural Network Schrödinger-Kohn) method. If `run_sk` is set to `True`, the model will be run using + the NNSK method. If + use_correction : Optional[str] + The `use_correction` parameter is an optional parameter that specifies whether to use correction + terms in the model. It can be set to either `None` or a string value. If it is set to `None`, the + model will not use any correction terms. If it is set to a string value + + Returns + ------- + The function `load_dptb_model` returns three variables: `allbonds`, `hamil_block`, and + `overlap_block`. + + ''' + ## create BaseStruct + structure_base =BaseStruct( + atom=ase.io.read(structure_tbtrans_file), + format='ase', + cutoff=apiHrk.apihost.model_config['bond_cutoff'], + proj_atom_anglr_m=apiHrk.apihost.model_config['proj_atom_anglr_m'], + proj_atom_neles=apiHrk.apihost.model_config['proj_atom_neles'], + onsitemode=apiHrk.apihost.model_config['onsitemode'], + time_symm=apiHrk.apihost.model_config['time_symm'] + ) + + apiHrk.update_struct(structure_base) + allbonds,hamil_block,overlap_block = apiHrk.get_HR() + + return allbonds,hamil_block,overlap_block + + def hamiltonian_get(self,allbonds:torch.tensor,hamil_block:torch.tensor,overlap_block:torch.tensor,Hamil_sisl,energy_unit_option:str): + '''The function `hamiltonian_get` takes in various parameters and calculates the Hamiltonian matrix + for a given set of bonds, storing the result in the `Hamil_sisl` matrix. + + Parameters + ---------- + allbonds + A torch.tensor containing information about the bonds in the system. Each element of the list + corresponds to a bond and contains information such as the indices of the atoms involved in the + bond and the displacement vector between the atoms. + hamil_block + The `hamil_block` parameter is a block of the Hamiltonian matrix. It is a tensor that contains the + values of the Hamiltonian matrix elements for a specific bond in the system. + overlap_block + The `overlap_block` parameter is a block matrix representing the overlap between orbitals in the + Hamiltonian. It is used to calculate the Hamiltonian matrix elements. + Hamil_sisl + Hamil_sisl is sisl.hamiltonian that represents the Hamiltonian matrix. The first two + dimensions correspond to the orbital indices, and the third dimension represents the supercell + indices. The Hamiltonian matrix elements are stored in this array. + Inertially, Hamil_sisl is a scipy.sparse.csr_matrix. + energy_unit_option + The `energy_unit_option` parameter is a string that specifies the unit of energy for the + calculation. It can be either "Hartree" or "eV". + + ''' + + if energy_unit_option=='Hartree': + unit_constant = 1 + + elif energy_unit_option=='eV': + unit_constant = 27.2107 + + else: + raise RuntimeError("energy_unit_option should be 'Hartree' or 'eV'") + + + # print(len(allbonds)) + # H_device.H[1000,1000]=1 + + x_max = abs(allbonds[:,-3].numpy()).max() + y_max = abs(allbonds[:,-2].numpy()).max() + z_max = abs(allbonds[:,-1].numpy()).max() + # print('x_max: ',x_max) + # print('y_max: ',y_max) + # print('z_max: ',z_max) + Hamil_sisl.set_nsc(a=2*abs(x_max)+1,b=2*abs(y_max)+1,c=2*abs(z_max)+1) + # set the number of super-cells in Hamiltonian object in sisl, which is based on allbonds results + + + for i in range(len(allbonds)): + # if i%100==0:print('bond_index: ',i) + orb_first_a = Hamil_sisl.geometry.a2o(allbonds[i,1]) + orb_last_a = Hamil_sisl.geometry.a2o(allbonds[i,1]+1) + orb_first_b = Hamil_sisl.geometry.a2o(allbonds[i,3]) + orb_last_b = Hamil_sisl.geometry.a2o(allbonds[i,3]+1) + + + if allbonds[i][-3:].equal(torch.tensor([0,0,0])): #allbonds[i,1] is atom index + + for orb_a in range(orb_first_a,orb_last_a): + for orb_b in range(orb_first_b,orb_last_b): + Hamil_sisl[orb_a,orb_b]=hamil_block[i].detach().numpy()[orb_a-orb_first_a,orb_b-orb_first_b]*unit_constant + # Hamil_sisl[orb_b,orb_a]=hamil_block[i].detach().numpy()[orb_b-orb_first_b,orb_a-orb_first_a]*unit_constant + Hamil_sisl[orb_b,orb_a]=np.conjugate(Hamil_sisl[orb_a,orb_b]) + else: + x = allbonds[i,-3].numpy().tolist() + y = allbonds[i,-2].numpy().tolist() + z = allbonds[i,-1].numpy().tolist() + # consistent with supercell setting:Set the number of super-cells in the `Lattice` object in sisl + # if abs(x) > 1 or abs(y) > 1 or abs(z) > 1: + + # print("Unexpected supercell index: ",[x,y,z]) + # print("Attention: the supercell setting may be too small to satisfy the nearest cell interaction, \ + # error in Lead self-energy calculation may occur.") + + for orb_a in range(orb_first_a,orb_last_a): + for orb_b in range(orb_first_b,orb_last_b): + H_value = hamil_block[i].detach().numpy()[orb_a-orb_first_a,orb_b-orb_first_b]*unit_constant + if H_value != 0: + Hamil_sisl[orb_a,orb_b,(x,y,z)]=H_value + Hamil_sisl[orb_b,orb_a,(-1*x,-1*y,-1*z)]=np.conjugate(Hamil_sisl[orb_a,orb_b,(x,y,z)]) + # Hamil_sisl[orb_b,orb_a,(-1*x,-1*y,-1*z)]=hamil_block[i].detach().numpy()[orb_b-orb_first_b,orb_a-orb_first_a]*unit_constant + + #TODO: At this stage, there is some problem using slice operation in sisl. I'm fixing it with the developer of sisl. + # I believe that the slice operation will be take soon. + + diff --git a/dptb/tests/data/test_negf/lead_R.nc b/dptb/tests/data/test_negf/lead_R.nc new file mode 100644 index 0000000000000000000000000000000000000000..aaffc5b710e2dd8864d7277d5c26402b2457f1ae GIT binary patch literal 61494 zcmeHQ3w%_?)t}8S5MV(-d04P4uZp~sM=2I@6G%dk5FSD0Q%hWug)AhyF%KT%S5{P7 z(9gDtl`6_-eH5*&#rLBh>*r7HSGB0ME!rp*srs|^-HKIFzH{!(f3lg~o0}zxB=!zW z=H5B;ICIXKb7sz+yLWS8UO}%OU+a;Toh=X!t1$I_nRAK){P^y?o%0Lx7G?SG%^Ggy zdHX2PieB>K6O|6_E$@3|2@hZOA0m8|cs&A*o2h>xm-C9=BAbhr7SCEx6RQr_)J4Lz z;#dMrpqNE*b-1QB5{u55G%gyhD=S|5E| zCry|nx=uX&%pUoR{UHprA;-$1%ozhYht579t=9W{ZIvs8kngEZ!b{ zXwW`!V3DkNx2lg3cC{bjST5=nUF8U873G$cSR-R1^2qY$%*!t>nKQq*v}9IpQNHk> zBM28S4Z$Kwaqhf)t9*_szqle&8=(K808wIAUO}KdLc?d`1eVp;1){OKM$qa|O(e9e zGCU!GFj4pkp)8B!eo1nRo>TH7RgrShhvb1-#uQ-WWcYR7852@z*j+)9wv2{C%$+AI zUZv`ThQ7{1!&5V(p_H*tcZ91`Y3Nln&}ob0u)DbA8?s`SY+TQ@o}Gn;lQW~Cl*P|0 zx~)m2VL;IUsk0nr_hagLJt>nEvAHPsgQj#z770~lzAy%?D zP1@IWcyh^RuIJAxoZGHCd-5a^`pe(b z(=csb&xammns*(@qwuJyOJ*w?y|C63ltYvKAMclUr^&?ndvNa;PYTiQObbjCKm5ff z|4L<=*7OP|uV0YBCEMLDJvhfZhDSyhacHY&PnK2Gpjx7>a=C_A94Gt~X>}z{{j6B5 zraTf2)zKO&x#M&TBUZeUnD_ON-7k@e%I#|{sf3_Kiwu2n(IgE;ToZ^b57gGvGO{dOSxJkf zXdoOas|eIbBX#3Bi#$){t%+r^)Z@d<@|3)yR%~4r z1gmRL;UDMAx*vJpd_-PpYFf0&BAL=ulm%#6D}$a1($)lm3X&tH^b}JL6SzxnITKJr zh#O~1N?*$YLm^Kln$pdzA@jLtMsqmjk6v)2Uj~WREO{O@mJ?5Nin& zK&~l~@UN+2NH6)|nlooEE-DaRetO!PY_X8qIbKXYG`No#MVI}>)ahgUi0xEqfLMNL z!wAaZ%bGiLr&tzyu-NqQ^WPTVrm_RY3kBYDMU)B*6z43Mw@&o7a{bp{zC}Dmx&1|c z&O`T!zgPtZ|7gbp;%w{kv8PV{jkudG2Z=YcNB>Un&TWAB(Y5Pe5(ll^B^UkIdm`7$ z9o(bpkhqL42a2-`FS>~?+|X(YSp4s*n~^Lg-pL@@p+#Ks0ZK%%kmi64lKW6pic`92 zDMxBaDHPAOmXR`|r?Zq3HQR*}1)e@uP_l(vDXhvLB9uJpRg9Eso>W%grRHF<_pP!J za{j$jwmwq7d({4Rqauj^>#W?gQ-m$KwQgX3vqa2{Rglrm897e*DaLp_Be^QuABc!; zdIa`omDRcsz^?^9jMm#99wf zHu|NC=9DZZWlRe}ioELPnUZd?$WwKdS`hN}ON)NCY4)u7MSKkdV_nlNX4@LY7Dk){ z^6yo%6I?U%G1JX_OnKy^eaw9`A42FD)rKuoIR8$|GN~S?8l2R@xj; z9^)M#vvof3v37Au<&i(graa3b2egm!GFOOgARsc_)TMw+0ha4)78nP4+nOPCogbcZdUydFNHtnNE^1 z@=lc50R{E9sSkv)eEhSMK1%Kp-hcGS{w}ZB&(2K>*U3bPvvAS&>FVd!otTVK!^yIb zM^#sC%N6$1By)w`9MN871ClWSP2dypzbr%#OD;}fPt>CPg8W7K#k2BDEl<=`A}f10 z%0qRb+HhT?OG+Pv64?xw9-@Hdse-ovX#{Q@D#zhm)iu0R4RD7@u)sawd*(_2?w`I- z6;uiC0YB@3rtr@%->*EJt;M9od*!Zz!n~ORhw-?hAM7h!+4yAK8;geFhVl4-k~dmf z3n!AthY96fyY51`R6IU7-+a+nY3(D*<<|IQkhAOIxbpbmdh!`e zx9u%!$pa{*QR^IuV*?iBNj1nvQ^&Efnis0A6@AkjW=QwILNlBGksi*k4D*?1qfCxu__1;*@+?egYjf`9=GFfA#3{Au0J(l|y{{S>1fu4PA0je(Bx= z71i(b%Q5p47pz-1{C_XrJ$y>Rhus?b>+;9%*u9j?oB1ogJ9;JI^XSF1H+}Zog@bkZ zvrgHug7IvcUzB5EljhF@Bkv)6hNr+^8?z}w{XRN*`^Hc5%L6$~6U-hKA@iO^oz07% zBCi(tc^2zxnu#Jnv1!Q3nL-&yGd9IXjt6)#2a6~rt#zN^$Qu~TX0EVAY8 z>#XktN`KC$J%;~Q(y$_Fzeoe58h1Y zH;ogFgTE3zw$+Y5858ayKZpM|hBLEr)n+l-ePh--GAS5qBL(g1%Hg5;$u*l!3V5r) ztX&-A*C@&=LNzPGZR@^aIh(fY`Qmg;R^C_Pzesul+Rm-RW8#1m&QY05Vs{#@WV-w8 zE7w|6ztQCgAFtd$r2d5$N1pC|R88z4+LZzWed@BhCmC3U}(L|sTGQd{bd<$=*8 zif3DgcfDCzH>xEHS7qRrLSD=7Hak6e#}g7v-V(f}N3Z-<)>lR8UA^*1S=f4KudX5- zsI3mw)P@7gV>MNwy0#@WKlpiO+P|hXt83fcpyCBOLdFU={o&ldEMAnK?;!@7YE6J&Em*EG39L#&$rDMvE zPTZpw=Cly9LJL-o`CKp~R}40uNO=Ru*ZyhgSC8~8+=9=EmWS6i`LPFDM(W?T;B#!5 zz0ZjLkuo0ezwo)AT*xToNq_34(6BMHW#{rxd|Z!0?#C-gn`YAFt9G`DkbdgFC0 zJa7Dq3wjr^r%Lw1>q{=s0$IrA8X3s~)0hKNYZoy(y&Ui%E#}pH`gORgPXQ6R28YLi$x=$Xbar5q$V$pER#!vJC|ew^B9}?1uD&Kb^85*J$Cgk(-9m(pC2EvcbVNw?D=5T$tI1>7^}E`6j^qert)_;i#b1r5wC_*96^znekW1nL)7f&A4idw(32=mrZ|Nn|FQy~ zKc|zYf*-Az^XHG0_a};bDZWV2X+C{IWj>|&8O1{so#u}_9+yAK`ns*%c#>gu zbq}UEgrbxEnMY-cDbA<3fZ~}H`TeQ0DK4RC&Yw)}k2@az9kj;vfXkoumACd@zqQsr zne)k9-#OWz*3KtWAGqVmZ2kRJ=TB?ptF850h;-*Nie(hbDTXO7r?`S*1;r|gQHn8& z)f7$sxZ`PG`{VNG%iV8hs;|xYaxu-3jTA4Tcqzrr6t_^koFci>Qsz%bv#0KOzFP8T z>SxoR*4FmqJ0#;vDV9<^pQ4lgicp!A6jxEKq-ffkTPg20ia(>clcGDG=KY;c`%RZW znaSIiz1}&MWIWmb<}`oI`I4#qG3Qe+>OT|zyOaKK$J2EE{O`se+w-FAY)%utsBfBy z$HV-r$PF$9Tne}pa4FzYz@@-fnF8G zItK+h2bS)9a{1IbD9|ymbmvRQ;DZwuE}xtz)G^>LpE?E~oUm~DKt&Fk1n4&2L(C@ zmhOCV`P4Zm&@wDtd2Sh=RPe+>rE$iRqQ_JXdv`e{sI@%oQ82!6)&@uSngoVo| zCkk~8xXY)G!3QTSIzOK>eP0+9nfhy7$&}yQ=it6SEa~9>bt&Ldz@>mo0ha7|DuDlIBMW3!FZ@ly-j1RczpXIDif9IzKEH%t${!Mi3GV zpBIa+s)*Ip){S2fDq9s;HGT#2#QI1}s4QF+j@H%s>7ya(K9kPB$;nMD8#8A7ME+t^ zW}o0;@^eT|!4CepH1EVq{Z;|@k4piU0&PQq-o6nH`|CbVTy${%2m20Rl9>C=uRr(a z4XYDqJv1@?;c<1(GZpdGl z(7!{IyZVwXy9xeQ!MR1Fhdv$bzOe5Nc}ei5?QaeC`+4BW>4$F8{TT3s34i?9z`@0% zuMR$V{=s7khkZI@QDx7u)BFEW_usVsEy&mSH-&FS`??&)*Zae(?s>6#&gu6a+vhv^ z3!epk0{y>n{o0Rm8ZHk`y{Tu<>HXsxZnTf_>Go0ospr=Ik?JQT{*Cp)m;4Lvxbt6k z-WS~ZKih{qeK@uo{84-fUbX$=hD&dFME47LeAm|5yHDL8-0<%H;_H5MAhF<$QQsXo z^cD0kW)Hrm8^II)$6nj>6M`2%cH+_zLtoMD$o97;N+Xvojc@sJ;-BkBR-GTZUDsFe zpaIwQ&nx{!#XN6dcm5iG!{cH_@WqRs+5hGTy-fTQcKF-0uls*JamL&WMxVcX=tU)k zA3gR&u;R(*J{h!d)#(_2U9d~!8~t`A?%4dd2`iu2mZ+}SxqZdkKhX7U#;5UD^@BD4 z_}jaY#`+qL`2xQ5KD_gaPl!+G-%Q_nFlW@DKI=7n^RIe~_YSB_?-Cq*S$RuNx5E!N!qt2Tf-jS1U-7L2bMDpk4S2$Y&zi8~?^Ex;Cjq{^S-bbQ zXMFMBy1oZFw722=+K|WH0PjTcUg}kZz6!?;N7gzj~d?NbQ%@kGI(BXD9pQWPfb* zFUBWS|F-t6h4!nJ_9?~wfximBt!Dcge<%CVTKlC0RnISDUJEDtHQJrIpj&TulEno6Oc#1o3_6#xaV6>4$Qgj zR^5*QPndAnFX&0kM~ts&{hb(>#vk?t>u9XkZ16+>x*gQtJL1Vl|5h}7H`bdr^)WtO zA9met{{p^LZ+c_H1rvV<{^`S#st7u^fkMaYY}e}G-JTmQ@Uw4|A76n+CQYM|J8iZd^umd(p_A8zwY0F2TeHq6?XeCC;Z_Lu-o5YzCv!y{^8F6 zU(EfPQ9lKL6F=NT-f?fHw{E)m{)2}a_h}~F;8Wa$8||mSA#Y|q+#Bs9Tc=KNC;Mg8 zkL!ljMTYntaeXKK>|~#8^tT=k?3eia;)!#@S707$`Gj2!B=j4_`EgCcANIA zY5z|4N!v$_pV5EN!T!O%KtJgH75q$@+7|qqN~-W{&l;! zXyf{5R;Os+!M`M+|FKSnJgbJZ{)L{xd;q+IK3~%M1@;Ge4)Y24H>HpEPsY>3wf_$G zNi?mmpD$@R`iGvxd_?=oKW6J+!oCvnigmK8FJRa0_6JV>J174e_RV(tGweI<_BS^E zKTRL(i|zJb3O`%_q)q#p{|au>SND%~1?DUCBg|Mvbt)o&dABQ^0rK2ppl+?y%n!OX|K;lD$+ zPL2NK8lWzW{!`$F{Y-%y?VHUDC;W~2PVGCvZTmO$qxpQ?PJin83;k)QzhS?yPB!gV z)Bd%7*Z8*f`4jN7vyVpqTJHb>><_UI!Fm_)ruC29{=g4pV}BBM&-ZQgr{<57{>Hu% z@{V;f#^jA8 zhE8$hBQ@!We58;E+?y%plgf{qaHsYS`;}7PXy0t!v@a?3js8>M8QQm$`bPV<_&NB0 z^!^8W1oEmHZ1wXY*axgLAx~)Ez|VX>B48h|&W60Fub=vS2lf&B2&`98ANWC!8~(Yt z4SrgG8F2X5G~8%Ef%%U9o%}ye{x>`yf`3iR$uD6#p02Q+D z)bOvFVByq0=98)zHyd!OZ`cR3x?}%Na{{ufe`_#1kl;=-Q=POM6m13V-tN(5E ztEPjKeQ>gGPWA(Qf&JU&^HlD5E8b%J=cNCg?30uIX`%hHwNI_p|Au`?d4Az!ADrx) zlJhv^-PCKaf8dLa{j;s_r2n1llN!J6`G^edSBibbxXl#$Pvz_Jp`1Y{CpfC<)X6>= a`ak3L4f~gIePjF@Z(Ot%KcjtH{Qe)3xf*By literal 0 HcmV?d00001 diff --git a/dptb/tests/data/test_tbtrans/best_nnsk_b3.600_c3.600_w0.300.json b/dptb/tests/data/test_tbtrans/best_nnsk_b3.600_c3.600_w0.300.json new file mode 100644 index 00000000..7fddb9d8 --- /dev/null +++ b/dptb/tests/data/test_tbtrans/best_nnsk_b3.600_c3.600_w0.300.json @@ -0,0 +1,122 @@ +{ + "onsite": { + "N-N-2s-2s-0": [ + 0.02462027035653591, + 0.007205560803413391 + ], + "N-N-2s-2p-0": [ + 0.008309782482683659, + -0.007032226305454969 + ], + "N-N-2p-2p-0": [ + 0.012606431730091572, + 0.010783562436699867 + ], + "N-N-2p-2p-1": [ + 0.0068643586710095406, + -0.011892829090356827 + ], + "N-B-2s-2s-0": [ + 0.041020166128873825, + -0.007834071293473244 + ], + "N-B-2s-2p-0": [ + 26.214815139770508, + -28.22139549255371 + ], + "N-B-2p-2p-0": [ + 0.2541739046573639, + 0.3701082468032837 + ], + "N-B-2p-2p-1": [ + -0.052932459861040115, + 0.03325718641281128 + ], + "B-N-2s-2s-0": [ + -0.10863762348890305, + -0.10621777176856995 + ], + "B-N-2s-2p-0": [ + 24.486974716186523, + 26.85447883605957 + ], + "B-N-2p-2p-0": [ + 0.13345032930374146, + -0.3127709925174713 + ], + "B-N-2p-2p-1": [ + -0.03844165802001953, + 0.00011800970969488844 + ], + "B-B-2s-2s-0": [ + -0.007172069512307644, + 0.007495054975152016 + ], + "B-B-2s-2p-0": [ + 0.004442985635250807, + 0.0030813836492598057 + ], + "B-B-2p-2p-0": [ + -0.0013882589992135763, + -6.591381679754704e-05 + ], + "B-B-2p-2p-1": [ + -0.009361814707517624, + -0.017272837460041046 + ] + }, + "hopping": { + "N-N-2s-2s-0": [ + 0.05967297405004501, + -0.21457917988300323 + ], + "N-N-2s-2p-0": [ + 0.042178086936473846, + 0.5767796635627747 + ], + "N-N-2p-2p-0": [ + 0.1008036881685257, + 0.5027011632919312 + ], + "N-N-2p-2p-1": [ + -0.005753065925091505, + -1.0040007829666138 + ], + "N-B-2s-2s-0": [ + 0.06605493277311325, + 2.485130786895752 + ], + "N-B-2s-2p-0": [ + -0.19711704552173615, + -1.884203314781189 + ], + "N-B-2p-2s-0": [ + -0.06678403168916702, + -2.8326284885406494 + ], + "N-B-2p-2p-0": [ + -0.2129121571779251, + 3.601222276687622 + ], + "N-B-2p-2p-1": [ + 0.034046269953250885, + -4.79115104675293 + ], + "B-B-2s-2s-0": [ + -0.061647042632102966, + -0.4071168601512909 + ], + "B-B-2s-2p-0": [ + 0.048091448843479156, + 0.8385105729103088 + ], + "B-B-2p-2p-0": [ + -0.11344866454601288, + 0.5754562020301819 + ], + "B-B-2p-2p-1": [ + 0.007624409161508083, + -0.7790966033935547 + ] + } +} \ No newline at end of file diff --git a/dptb/tests/data/test_tbtrans/negf_tbtrans.json b/dptb/tests/data/test_tbtrans/negf_tbtrans.json new file mode 100644 index 00000000..c4d6b748 --- /dev/null +++ b/dptb/tests/data/test_tbtrans/negf_tbtrans.json @@ -0,0 +1,90 @@ +{ + "common_options": { + "onsitemode": "strain", + "onsite_cutoff": 1.6, + "bond_cutoff": 3.6, + "env_cutoff": 3.5, + "atomtype": [ + "N", + "B" + ], + "proj_atom_neles": { + "N": 5, + "B": 3 + }, + "proj_atom_anglr_m": { + "N": [ + "2s", + "2p" + ], + "B": [ + "2s", + "2p" + ] + } + }, + "model_options": { + "sknetwork": { + "sk_hop_nhidden": 1, + "sk_onsite_nhidden": 1 + }, + "skfunction": { + "sk_cutoff": 3.6, + "sk_decay_w": 0.3 + } + }, + "structure":"./test_hBN_zigzag_struct.xyz", + "task_options": + { + "task": "tbtrans_negf", + "scf": true, + "block_tridiagonal": false, + "ele_T": 500, + "unit": "Hartree", + "scf_options":{ + "mode": "PDIIS", + "mixing_period": 3, + "step_size": 0.05, + "n_history": 6, + "abs_err": 1e-6, + "rel_err": 1e-4, + "max_iter": 100 + }, + "stru_options":{ + "pbc":[false, true, false], + "kmesh":[1,1,1], + "device":{ + "id":"8-12", + "sort": true + }, + "lead_L":{ + "id":"0-8", + "voltage":0.0 + }, + "lead_R":{ + "id":"12-20", + "voltage":0.0 + } + }, + "poisson_options": { + "solver": "fmm", + "err": 1e-5 + }, + "sgf_solver": "Sancho-Rubio", + "espacing": 0.1, + "emin": -2, + "emax": 2, + "e_fermi": -9.874357223510742, + "density_options":{ + "method": "Ozaki" + }, + "eta_lead":1e-5, + "eta_device":0.0, + "out_dos": true, + "out_tc": true, + "out_ldos": true, + "out_current_nscf": true, + "out_density": true, + "out_lcurrent": true + } +} diff --git a/dptb/tests/data/test_tbtrans/test_hBN_zigzag_struct.xyz b/dptb/tests/data/test_tbtrans/test_hBN_zigzag_struct.xyz new file mode 100644 index 00000000..f05a9679 --- /dev/null +++ b/dptb/tests/data/test_tbtrans/test_hBN_zigzag_struct.xyz @@ -0,0 +1,22 @@ +20 +Lattice="30.0 0.0 0.0 0.0 4.337055133 0.0 0.0 0.0 12.519999742500001" Properties=species:S:1:pos:R:3 pbc="T T T" +N 15.00000000 3.97563387 0.62599999 +B 15.00000000 2.52994883 0.62599999 +N 15.00000000 1.80710631 1.87799996 +B 15.00000000 0.36142126 1.87799996 +N 15.00000000 3.97563387 3.12999994 +B 15.00000000 2.52994883 3.12999994 +N 15.00000000 1.80710631 4.38199991 +B 15.00000000 0.36142126 4.38199991 +N 15.00000000 3.97563387 5.63399988 +B 15.00000000 2.52994883 5.63399988 +N 15.00000000 1.80710631 6.88599986 +B 15.00000000 0.36142126 6.88599986 +N 15.00000000 3.97563387 8.13799983 +B 15.00000000 2.52994883 8.13799983 +N 15.00000000 1.80710631 9.38999981 +B 15.00000000 0.36142126 9.38999981 +N 15.00000000 3.97563387 10.64199978 +B 15.00000000 2.52994883 10.64199978 +N 15.00000000 1.80710631 11.89399976 +B 15.00000000 0.36142126 11.89399976 diff --git a/dptb/tests/test_tbtrans_init.py b/dptb/tests/test_tbtrans_init.py new file mode 100644 index 00000000..7bb49009 --- /dev/null +++ b/dptb/tests/test_tbtrans_init.py @@ -0,0 +1,106 @@ +import pytest +from dptb.nnops.apihost import NNSKHost +from dptb.plugins.init_nnsk import InitSKModel +from dptb.nnops.NN2HRK import NN2HRK +from dptb.postprocess.tbtrans_init import TBTransInputSet +from dptb.utils.tools import j_loader +import numpy as np + +@pytest.fixture(scope='session', autouse=True) +def root_directory(request): + return str(request.config.rootdir) + + +def test_tbtrans_init(root_directory): + + # check whether sisl is installed: if not, skip this test + try: + import sisl + except: + pytest.skip('sisl is not installed which is necessary for TBtrans Input Generation. Therefore, skipping test_tbtrans_init.') + + model_ckpt = f'{root_directory}/dptb/tests/data/test_tbtrans/best_nnsk_b3.600_c3.600_w0.300.json' + config = f'{root_directory}/dptb/tests/data/test_tbtrans/negf_tbtrans.json' + apihost = NNSKHost(checkpoint=model_ckpt, config=config) + apihost.register_plugin(InitSKModel()) + apihost.build() + apiHrk = NN2HRK(apihost=apihost, mode='nnsk') + + run_opt = { + "run_sk": True, + "init_model":model_ckpt, + "results_path":f'{root_directory}/dptb/tests/data/test_tbtrans/', + "structure":f'{root_directory}/dptb/tests/data/test_tbtrans/test_hBN_zigzag_struct.xyz', + "log_path": '/data/DeepTB/dptb_Zjj/DeePTB/dptb/tests/data/test_tbtrans/output', + "log_level": 5, + "use_correction":False + } + + jdata = j_loader(config) + jdata = jdata['task_options'] + tbtrans_hBN = TBTransInputSet(apiHrk=apiHrk,run_opt=run_opt, jdata=jdata) + + # check _orbitals_name_get + element_orbital_name = tbtrans_hBN._orbitals_name_get(['2s', '2p']) + assert element_orbital_name == ['2s', '2py', '2pz', '2px'] + element_orbital_name = tbtrans_hBN._orbitals_name_get(['3s', '3p', 'd*']) + assert element_orbital_name == ['3s', '3py', '3pz', '3px', 'dxy*', 'dyz*', 'dz2*', 'dxz*', 'dx2-y2*'] + + + + tbtrans_hBN.hamil_get_write(write_nc=False) + assert (tbtrans_hBN.allbonds_all[0].detach().numpy()-np.array([5, 0, 5, 0, 0, 0, 0])).max()<1e-5 + assert (tbtrans_hBN.allbonds_all[50].detach().numpy()-np.array([ 5, 2, 5, 18, 0, 0, -1])).max()<1e-5 + assert (tbtrans_hBN.hamil_block_all[0].detach().numpy()- np.array([[-0.73303634, 0. , 0. , 0. ], + [ 0. , 0.04233637, 0. , 0. ], + [ 0. , 0. , 0.04233636, 0. ], + [ 0. , 0. , 0. , -0.27156556]])).max()<1e-5 + assert (tbtrans_hBN.hamil_block_all[50].detach().numpy()-np.array([[-0.03164609, -0. , 0.02028139, -0. ], + [ 0. , 0.00330366, 0. , 0. ], + [-0.02028139, 0. , -0.05393751, 0. ], + [ 0. , 0. , 0. , 0.00330366]])).max()<1e-5 + assert (tbtrans_hBN.allbonds_lead_L[0].detach().numpy()-np.array([5, 0, 5, 0, 0, 0, 0])).max()<1e-5 + assert (tbtrans_hBN.allbonds_lead_L[50].detach().numpy()-np.array([5, 4, 7, 7, 0, 0, 0])).max()<1e-5 + assert (tbtrans_hBN.hamil_block_lead_L[0].detach().numpy()-np.array([[-0.73303634, 0. , 0. , 0. ], + [ 0. , 0.04233637, 0. , 0. ], + [ 0. , 0. , 0.04233636, 0. ], + [ 0. , 0. , 0. , -0.27156556]])).max()<1e-5 + assert (tbtrans_hBN.hamil_block_lead_L[50].detach().numpy()-np.array([[ 0.1145315 , -0.06116847, 0.10594689, 0. ], + [ 0.15539739, -0.04634972, 0.22751862, 0. ], + [-0.26915616, 0.22751862, -0.30906558, 0. ], + [-0. , 0. , 0. , 0.08500822]])).max()<1e-5 + + + # tbtrans_hBN.hamil_write() + # check the hamiltonian through Hk at Gamma and M + H_lead = tbtrans_hBN.H_lead_L + G_eigs = sisl.BandStructure(H_lead, [[0., 0.,0.]], 1, ["G"]) + M_eigs = sisl.BandStructure(H_lead, [[0, 0.5 ,0 ]], 1, ["M"]) + Ef = -9.874358177185059 + G_eigs = G_eigs.apply.array.eigh() -Ef + M_eigs = M_eigs.apply.array.eigh() -Ef + + G_eigs_right = np.array([[-19.95362763, -17.10774579, -17.10774579, -16.9118761 , + -11.20609829, -10.01050689, -10.01050689, -8.11443067, + -8.11443067, -7.60442045, -6.6550603 , -3.79211664, + -3.79211663, -3.22863532, -3.22863532, -3.17535758, + 3.18703263, 3.24031037, 3.24031037, 6.45306332, + 7.70583557, 7.91107113, 10.91058699, 10.91058699, + 23.64785516, 23.64785516, 28.30755414, 28.30755428, + 28.65719263, 30.78452851, 33.25399887, 33.25399887]]) + + M_eigs_right = np.array([[-18.69653568, -18.69653568, -16.91187582, -16.91187582, + -11.20609828, -11.20609828, -7.44726991, -7.44726991, + -6.65506047, -6.65506047, -5.79252308, -5.79252308, + -5.2193769 , -5.2193769 , -3.17535758, -3.17535758, + 3.18703263, 3.18703263, 5.84906816, 5.84906816, + 7.74726616, 7.74726616, 7.91107121, 7.91107121, + 28.12912079, 28.12912079, 28.65719227, 28.65719227, + 29.54182975, 29.54182975, 30.78452877, 30.78452877]]) + + assert (G_eigs[0]-G_eigs_right[0]).max()<1e-5 + assert (M_eigs[0]-M_eigs_right[0]).max()<1e-5 + + + + diff --git a/dptb/utils/argcheck.py b/dptb/utils/argcheck.py index 557acbdb..8b4e19b7 100644 --- a/dptb/utils/argcheck.py +++ b/dptb/utils/argcheck.py @@ -430,10 +430,66 @@ def task_options(): Argument("FS3D", dict, FS3D()), Argument("write_sk", dict, write_sk()), Argument("ifermi", dict, ifermi()), - Argument("negf", dict, negf()) + Argument("negf", dict, negf()), + Argument("tbtrans_negf", dict, tbtrans_negf()) ],optional=True, default_tag="band", doc=doc_task) +def tbtrans_negf(): + doc_scf = "" + doc_block_tridiagonal = "" + doc_ele_T = "" + doc_unit = "" + doc_scf_options = "" + doc_stru_options = "" + doc_poisson_options = "" + doc_sgf_solver = "" + doc_espacing = "" + doc_emin = "" + doc_emax = "" + doc_e_fermi = "" + doc_eta_lead = "" + doc_eta_device = "" + doc_out_dos = "" + doc_out_tc = "" + doc_out_current = "" + doc_out_current_nscf = "" + doc_out_ldos = "" + doc_out_density = "" + doc_out_lcurrent = "" + doc_density_options = "" + doc_out_potential = "" + + return [ + Argument("scf", bool, optional=True, default=False, doc=doc_scf), + Argument("block_tridiagonal", bool, optional=True, default=False, doc=doc_block_tridiagonal), + Argument("ele_T", [float, int], optional=False, doc=doc_ele_T), + Argument("unit", str, optional=True, default="Hartree", doc=doc_unit), + Argument("scf_options", dict, optional=True, default={}, sub_fields=[], sub_variants=[scf_options()], doc=doc_scf_options), + Argument("stru_options", dict, optional=False, sub_fields=stru_options(), doc=doc_stru_options), + Argument("poisson_options", dict, optional=True, default={}, sub_fields=[], sub_variants=[poisson_options()], doc=doc_poisson_options), + Argument("sgf_solver", str, optional=True, default="Sancho-Rubio", doc=doc_sgf_solver), + Argument("espacing", [int, float], optional=False, doc=doc_espacing), + Argument("emin", [int, float], optional=False, doc=doc_emin), + Argument("emax", [int, float], optional=False, doc=doc_emax), + Argument("e_fermi", [int, float], optional=False, doc=doc_e_fermi), + Argument("density_options", dict, optional=True, default={}, sub_fields=[], sub_variants=[density_options()], doc=doc_density_options), + Argument("eta_lead", [int, float], optional=True, default=1e-5, doc=doc_eta_lead), + Argument("eta_device", [int, float], optional=True, default=0., doc=doc_eta_device), + Argument("out_dos", bool, optional=True, default=False, doc=doc_out_dos), + Argument("out_tc", bool, optional=True, default=False, doc=doc_out_tc), + Argument("out_density", bool, optional=True, default=False, doc=doc_out_density), + Argument("out_potential", bool, optional=True, default=False, doc=doc_out_potential), + Argument("out_current", bool, optional=True, default=False, doc=doc_out_current), + Argument("out_current_nscf", bool, optional=True, default=False, doc=doc_out_current_nscf), + Argument("out_ldos", bool, optional=True, default=False, doc=doc_out_ldos), + Argument("out_lcurrent", bool, optional=True, default=False, doc=doc_out_lcurrent) + ] + + + + + def negf(): doc_scf = "" doc_block_tridiagonal = ""