Source code for triqs_dft_tools.converters.plovasp.plotools


################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Ferrero, O. Parcollet
#
# DFT tools: Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# PLOVasp: Copyright (C) 2015 by O. E. Peil
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
r"""
    plovasp.plotools
    ================

    Set of routines for processing and outputting PLOs.

    This is the main module containing routines responsible for checking
    the consistency of the input data, generation of projected localized
    orbitals (PLOs) out of raw VASP projectors, and outputting data
    required by DFTTools.

    The first step of PLO processing is to select subsets of projectors
    corresponding to PLO groups. Each group contains a set of shells. Each
    projector shell is represented by an object 'ProjectorShell' that contains
    an array of projectors and information on the shell itself (orbital number,
    ions, etc.). 'ProjectorShell's are contained in both a list of shells
    (according to the original list as read from config-file) and in a
    'ProjectorGroup' object, the latter also providing information about the
    energy window.

    Order of operations:
     - transform projectors (all bands) in each shell
     - select transformed shell projectors for a given group within the window
     - orthogonalize if necessary projectors within a group by performing
       the following operations for each k-point:
       * combine all projector shells into a single array
       * orthogonalize the array
       * distribute back the arrays assuming that the order is preserved

"""
import itertools as it
import numpy as np
from .proj_group import ProjectorGroup
from .proj_shell import ProjectorShell
from .proj_shell import ComplementShell

np.set_printoptions(suppress=True)

# 'simplejson' is supposed to be faster than 'json' in stdlib.
try:
    import simplejson as json
except ImportError:
    import json

[docs] def issue_warning(message): """ Issues a warning. """ print() print(" !!! WARNING !!!: " + message) print()
################################################################################ # check_data_consistency() ################################################################################
[docs] def check_data_consistency(pars, el_struct): """ Check the consistency of the VASP data. """ # Check that ions inside each shell are of the same sort for sh in pars.shells: max_ion_index = max([max(gr) for gr in sh['ions']['ion_list']]) assert max_ion_index <= el_struct.natom, "Site index in the projected shell exceeds the number of ions in the structure" ion_list = list(it.chain(*sh['ions']['ion_list'])) sorts = set([el_struct.type_of_ion[io] for io in ion_list]) assert len(sorts) == 1, "Each projected shell must contain only ions of the same sort" # Check that ion and orbital lists in shells match those of projectors lshell = sh['lshell'] for ion in ion_list: for par in el_struct.proj_params: if par['isite'] - 1 == ion and par['l'] == lshell: break else: errmsg = "Projector for isite = %s, l = %s does not match PROJCAR"%(ion + 1, lshell) raise Exception(errmsg)
################################################################################ # # generate_plo() # ################################################################################
[docs] def generate_plo(conf_pars, el_struct): """ Parameters ---------- conf_pars (dict) : dictionary of input parameters (from conf-file) el_struct : ElectronicStructure object """ check_data_consistency(conf_pars, el_struct) proj_raw = el_struct.proj_raw try: efermi = conf_pars.general['efermi'] except (KeyError, AttributeError): efermi = el_struct.efermi # eigvals(nktot, nband, ispin) are defined with respect to the Fermi level eigvals = el_struct.eigvals - efermi # check if at least one shell is correlated assert np.any([shell['corr'] for shell in conf_pars.shells]), 'at least one shell has be CORR = True' nshell = len(conf_pars.shells) print() print(" Generating %i shell%s..."%(nshell, '' if nshell == 1 else 's')) pshells = [] for sh_par in conf_pars.shells: pshell = ProjectorShell(sh_par, proj_raw, el_struct.proj_params, el_struct.kmesh, el_struct.structure, el_struct.nc_flag) print() print(" Shell : %s"%(pshell.user_index)) print(" Orbital l : %i"%(pshell.lorb)) print(" Number of ions: %i"%(pshell.nion)) print(" Dimension : %i"%(pshell.ndim)) print(" Correlated : %r"%(pshell.corr)) print(" Ion sort : %r"%(pshell.ion_sort)) pshells.append(pshell) pgroups = [] for gr_par in conf_pars.groups: pgroup = ProjectorGroup(gr_par, pshells, eigvals) pgroup.orthogonalize() if pgroup.complement: pgroup.calc_complement(eigvals) if conf_pars.general['hk']: pgroup.calc_hk(eigvals) #testout = 'hk.out.h5' #from h5 import HDFArchive #with HDFArchive(testout, 'w') as h5test: # h5test['hk'] = pgroup.hk # DEBUG output print("Density matrix:") nimp = 0.0 ov_all = [] for ish in pgroup.ishells: if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): print(" Shell %i"%(ish + 1)) dm_all, ov_all_ = pshells[ish].density_matrix(el_struct) ov_all.append(ov_all_[0]) spin_fac = 2 if dm_all.shape[0] == 1 else 1 for io in range(dm_all.shape[1]): print(" Site %i"%(io + 1)) dm = spin_fac * dm_all[:, io, : ,:].sum(0) for row in dm: print(''.join(map("{0:14.7f}".format, row))) ndm = dm.trace() if pshells[ish].corr: nimp += ndm print(" trace: ", ndm) print() print(" Impurity density:", nimp) print() print("Overlap:") for io, ov in enumerate(ov_all): print(" Site %i"%(io + 1)) print(ov[0,...]) print() print("Local Hamiltonian:") for ish in pgroup.ishells: if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): print(" Shell %i"%(ish + 1)) loc_ham = pshells[pgroup.ishells[ish]].local_hamiltonian(el_struct) for io in range(loc_ham.shape[1]): print(" Site %i (real | complex part)"%(io + 1)) for row in loc_ham[:, io, :, :].sum(0): print(''.join(map("{0:14.7f}".format, row.real))+' |'+''.join(map("{0:14.7f}".format, row.imag))) # END DEBUG output if 'dosmesh' in conf_pars.general: print() print("Evaluating DOS...") mesh_pars = conf_pars.general['dosmesh'] if np.isnan(mesh_pars['emin']): dos_emin = pgroup.emin dos_emax = pgroup.emax else: dos_emin = mesh_pars['emin'] dos_emax = mesh_pars['emax'] n_points = mesh_pars['n_points'] emesh = np.linspace(dos_emin, dos_emax, n_points) for ish in pgroup.ishells: if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell) or True: print(" Shell %i"%(ish + 1)) dos = pshells[pgroup.ishells[ish]].density_of_states(el_struct, emesh) de = emesh[1] - emesh[0] ntot = (dos[1:,...] + dos[:-1,...]).sum(0) / 2 * de print(" Total number of states:", ntot) for io in range(dos.shape[2]): np.savetxt('pdos_%i_%i.dat'%(ish,io), np.vstack((emesh.T, dos[:, 0, io, :].T)).T) pgroups.append(pgroup) return pshells, pgroups
################################################################################ # # output_as_text # ################################################################################
[docs] def output_as_text(pars, el_struct, pshells, pgroups): """ Output all information necessary for the converter as text files. """ ctrl_output(pars, el_struct, len(pgroups)) plo_output(pars, el_struct, pshells, pgroups) if pars.general['hk']: hk_output(pars, el_struct, pgroups)
# TODO: k-points with weights should be stored once and for all ################################################################################ # # kpoints_output # ################################################################################
[docs] def kpoints_output(basename, el_struct): """ Outputs k-point data into a text file. """ kmesh = el_struct.kmesh fname = basename + '.kpoints' with open(fname, 'wt') as f: f.write("# Number of k-points: nktot\n") nktot = kmesh['nktot'] f.write("%i\n"%(nktot)) # TODO: add the output of reciprocal lattice vectors f.write("# List of k-points with weights\n") for ik in range(nktot): kx, ky, kz = kmesh['kpoints'][ik, :] kwght = kmesh['kweights'][ik] f.write("%15.10f%15.10f%15.10f%20.10f\n"%(kx, ky, kz, kwght)) # Check if there are tetrahedra defined and if they are, output them try: ntet = kmesh['ntet'] volt = kmesh['volt'] f.write("\n# Number of tetrahedra and volume: ntet, volt\n") f.write("%i %s\n"%(ntet, volt)) f.write("# List of tetrahedra: imult, ik1, ..., ik4\n") for it in range(ntet): f.write(' '.join(map("{0:d}".format, *kmesh['itet'][it, :])) + '\n') except KeyError: pass
################################################################################ # # ctrl_output # ################################################################################
[docs] def ctrl_output(conf_pars, el_struct, ng): """ Outputs a ctrl-file. Control file format """""""""""""""""""""""""""""" Filename '<namebase>.ctrl'. Contains the data shared between all shells. The JSON-header consists of the following elements: * *nk*: number of `k`-points * *ns*: number of spin channels * *nc_flag*: collinear/noncollinear case (False/True) * *ng*: number of projector groups * Symmetry information (list of symmetry operations) * *efermi*: Fermi level (optional) * Lattice information """ ctrl_fname = conf_pars.general['basename'] + '.ctrl' head_dict = {} # TODO: Add output of tetrahedra # Construct the header dictionary head_dict['ngroups'] = ng head_dict['nk'] = el_struct.kmesh['nktot'] head_dict['ns'] = el_struct.nspin head_dict['kvec1'] = list(el_struct.structure['kpt_basis'][:,0]) head_dict['kvec2'] = list(el_struct.structure['kpt_basis'][:,1]) head_dict['kvec3'] = list(el_struct.structure['kpt_basis'][:,2]) head_dict['nc_flag'] = 1 if el_struct.nc_flag else 0 # head_dict['efermi'] = conf_pars.general['efermi'] # We probably don't need Efermi header = json.dumps(head_dict, indent=4, separators=(',', ': ')) print(" Storing ctrl-file...") with open(ctrl_fname, 'wt') as f: f.write(header + "\n") f.write("#END OF HEADER\n") f.write("# k-points and weights\n") labels = ['kx', 'ky', 'kz', 'kweight'] out = "".join([s.center(15) for s in labels]) f.write("#" + out + "\n") for ik, kp in enumerate(el_struct.kmesh['kpoints']): tmp1 = "".join(map("{0:15.10f}".format, kp)) out = tmp1 + "{0:16.10f}".format(el_struct.kmesh['kweights'][ik]) f.write(out + "\n") f.write("# k-points and weights cartesian\n") labels = ['kx', 'ky', 'kz'] out = "".join([s.center(15) for s in labels]) f.write("#" + out + "\n") for ik, kp in enumerate(el_struct.kmesh['kpoints_cart']): out = "".join(map("{0:15.10f}".format, kp)) f.write(out + "\n")
################################################################################ # # plo_output # ################################################################################
[docs] def plo_output(conf_pars, el_struct, pshells, pgroups): """ Outputs PLO groups into text files. Filenames are defined by <basename> that is passed from config-file. All necessary general parameters are stored in a file '<basename>.ctrl'. Each group is stored in a '<basename>.plog<Ng>' file. The format is the following: | # Energy window: emin, emax | ib_min, ib_max | nelect | # Eigenvalues | isp, ik1, kx, ky, kz, kweight | ib1, ib2 | eig1 | eig2 | ... | eigN | ik2, kx, ky, kz, kweight | ... | # Projected shells | Nshells | # Shells: <shell indices> | # Shell <1> | Shell 1 | ndim | # complex arrays: plo(ns, nion, ndim, nb) | ... | # Shells: <shell indices> | # Shell <2> | Shell 2 | ... """ for ig, pgroup in enumerate(pgroups): plo_fname = conf_pars.general['basename'] + '.pg%i'%(ig + 1) print(" Storing PLO-group file '%s'..."%(plo_fname)) head_dict = {} head_dict['nb_max'] = pgroup.nb_max if 'bands' in conf_pars.groups[ig]: head_dict['bandwindow'] = (pgroup.ib_min, pgroup.ib_max) else: head_dict['ewindow'] = (pgroup.emin, pgroup.emax) # Number of electrons within the window head_dict['nelect'] = pgroup.nelect_window(el_struct) print(" Density within window:", head_dict['nelect']) head_shells = [] for ish in pgroup.ishells: shell = pgroup.shells[ish] sh_dict = {} sh_dict['shell_index'] = ish sh_dict['lorb'] = shell.lorb sh_dict['ndim'] = shell.ndim sh_dict['corr'] = shell.corr # Convert ion indices from the internal representation (starting from 0) # to conventional VASP representation (starting from 1) ion_output = [io + 1 for io in shell.ion_list] # Derive sorts from equivalence classes sh_dict['ion_list'] = ion_output sh_dict['ion_sort'] = shell.ion_sort # TODO: add the output of transformation matrices head_shells.append(sh_dict) head_dict['shells'] = head_shells header = json.dumps(head_dict, indent=4, separators=(',', ': ')) with open(plo_fname, 'wt') as f: f.write(header + "\n") f.write("#END OF HEADER\n") # Eigenvalues within the window if 'bands' in conf_pars.groups[ig]: f.write("# Eigenvalues within the band window: %s, %s\n"%(pgroup.ib_min+1, pgroup.ib_max+1)) else: f.write("# Eigenvalues within the energy window: %s, %s\n"%(pgroup.emin, pgroup.emax)) nk, nband, ns_band = el_struct.eigvals.shape for isp in range(ns_band): f.write("# is = %i\n"%(isp + 1)) for ik in range(nk): ib1, ib2 = pgroup.ib_win[ik, isp, 0], pgroup.ib_win[ik, isp, 1] # Output band indices in Fortran convention! f.write(" %i %i\n"%(ib1 + 1, ib2 + 1)) for ib in range(ib1, ib2 + 1): eigv_ef = el_struct.eigvals[ik, ib, isp] - el_struct.efermi f_weight = el_struct.ferw[isp, ik, ib] f.write("%13.8f %12.7f\n"%(eigv_ef, f_weight)) # Projected shells f.write("# Projected shells\n") f.write("# Shells: %s\n"%(pgroup.ishells)) for ish in pgroup.ishells: shell = pgroup.shells[ish] f.write("# Shell %i\n"%(ish)) nion, ns, nk, nlm, nb = shell.proj_win.shape for isp in range(ns): f.write("# is = %i\n"%(isp + 1)) for ik in range(nk): f.write("# ik = %i\n"%(ik + 1)) for ion in range(nion): for ilm in range(nlm): ib1, ib2 = pgroup.ib_win[ik, isp, 0], pgroup.ib_win[ik, isp, 1] ib_win = ib2 - ib1 + 1 for ib in range(ib_win): p = shell.proj_win[ion, isp, ik, ilm, ib] f.write("{0:16.10f}{1:16.10f}\n".format(p.real, p.imag)) f.write("\n")
################################################################################ # # plo_output # ################################################################################
[docs] def hk_output(conf_pars, el_struct, pgroups): """ Outputs HK into text file. Filename is defined by <basename> that is passed from config-file. The Hk for each groups is stored in a '<basename>.hk<Ng>' file. The format is similar as defined in the Hk dft_tools format, but does not store info about correlated shells and irreps nk # number of k-points n_el # electron density n_sh # number of total atomic shells at sort l dim # atom, sort, l, dim at sort l dim # atom, sort, l, dim After these header lines, the file has to contain the Hamiltonian matrix in orbital space. The standard convention is that you give for each k-point first the matrix of the real part, then the matrix of the imaginary part, and then move on to the next k-point. """ for ig, pgroup in enumerate(pgroups): hk_fname = conf_pars.general['basename'] + '.hk%i'%(ig + 1) print(" Storing HK-group file '%s'..."%(hk_fname)) head_shells = [] for ish in pgroup.ishells: shell = pgroup.shells[ish] ion_output = [io + 1 for io in shell.ion_list] for iion in ion_output: sh_dict = {} sh_dict['shell_index'] = ish sh_dict['lorb'] = shell.lorb sh_dict['ndim'] = shell.ndim # Convert ion indices from the internal representation (starting from 0) # to conventional VASP representation (starting from 1) # Derive sorts from equivalence classes sh_dict['ion_list'] = ion_output sh_dict['ion_sort'] = shell.ion_sort head_shells.append(sh_dict) with open(hk_fname, 'wt') as f: # Eigenvalues within the window nk, nband, ns_band = el_struct.eigvals.shape f.write('%i # number of kpoints\n'%nk) f.write('{0:0.4f} # electron density\n'.format(pgroup.nelect_window(el_struct))) f.write('%i # number of shells\n'%len(head_shells)) for head in head_shells: f.write('%i %i %i %i # atom sort l dim\n'%(head['ion_list'][0],head['ion_sort'][0],head['lorb'],head['ndim'])) norbs = pgroup.hk.shape[2] for isp in range(ns_band): for ik in range(nk): for io in range(norbs): for iop in range(norbs): f.write(" {0:14.10f}".format(pgroup.hk[isp,ik,io,iop].real)) f.write("\n") for io in range(norbs): for iop in range(norbs): f.write(" {0:14.10f}".format(pgroup.hk[isp,ik,io,iop].imag)) f.write("\n")