mirror of
https://github.com/triqs/dft_tools
synced 2025-05-07 07:34:58 +02:00
new VASP 6.5.0 hdf5 support capability (#262)
* [vasp] read fermi from vasp h5 archive and write deltaN to vasp h5 archive * Reading vasptriqs.h5 file; reading GAMMA file for non-collinear full csc * fix read error for tetrahedron data and eigs properties * write deltaN only to vasp h5 if present, write to GAMMA text file only for old interface * support symmetries at DFT level in VASP * use IBZ kpoint list for writing GAMMA vaspgamma.h5 * make vaspgamma.h5 default * write n_k_ibz into dft_misc_input so in can be used inside of sumk * add new svo converter test --------- Co-authored-by: Dario Fiore Mosca <fiosca@Darios-MacBook-Air.local>
This commit is contained in:
parent
54918d74ba
commit
2a8acf45a3
@ -39,6 +39,7 @@ class ElectronicStructure:
|
||||
|
||||
- *natom* (int) : total number of atoms
|
||||
- *nktot* (int) : total number of `k`-points
|
||||
- *nkibz* (int) : number of `k`-points in IBZ
|
||||
- *nband* (int) : total number of bands
|
||||
- *nspin* (int) : spin-polarization
|
||||
- *nc_flag* (True/False) : non-collinearity flag
|
||||
@ -58,8 +59,9 @@ class ElectronicStructure:
|
||||
self.natom = vasp_data.poscar.nq
|
||||
self.type_of_ion = vasp_data.poscar.type_of_ion
|
||||
self.nktot = vasp_data.kpoints.nktot
|
||||
self.nkibz = vasp_data.kpoints.nkibz
|
||||
|
||||
self.kmesh = {'nktot': self.nktot}
|
||||
self.kmesh = {'nktot': self.nktot, 'nkibz': self.nkibz}
|
||||
self.kmesh['kpoints'] = vasp_data.kpoints.kpts
|
||||
# VASP.6.
|
||||
self.nc_flag = vasp_data.plocar.nc_flag
|
||||
@ -86,21 +88,15 @@ class ElectronicStructure:
|
||||
|
||||
# Check that the number of k-points is the same in all files
|
||||
_, ns_plo, nk_plo, nb_plo = vasp_data.plocar.plo.shape
|
||||
assert nk_plo == self.nktot, "PLOCAR is inconsistent with IBZKPT (number of k-points)"
|
||||
assert nk_plo == self.nktot, "PLOCAR is inconsistent with IBZKPT (number of k-points). If you run VASP with symmetry make sure to use the h5 interface of the converter, i.e. have the locproj information written to vaspout.h5"
|
||||
|
||||
# FIXME: Reading from EIGENVAL is obsolete and should be
|
||||
# removed completely.
|
||||
# if not vasp_data.eigenval.eigs is None:
|
||||
if False:
|
||||
if vasp_data.eigenval.eigs is not None:
|
||||
print("eigvals from EIGENVAL")
|
||||
self.eigvals = vasp_data.eigenval.eigs
|
||||
self.ferw = vasp_data.eigenval.ferw.transpose((2, 0, 1))
|
||||
|
||||
nk_eig = vasp_data.eigenval.nktot
|
||||
assert nk_eig == self.nktot, "PLOCAR is inconsistent with EIGENVAL (number of k-points)"
|
||||
|
||||
# Check that the number of band is the same in PROJCAR and EIGENVAL
|
||||
assert nb_plo == self.nband, "PLOCAR is inconsistent with EIGENVAL (number of bands)"
|
||||
self.efermi = vasp_data.doscar.efermi
|
||||
else:
|
||||
print("eigvals from LOCPROJ")
|
||||
self.eigvals = vasp_data.plocar.eigs
|
||||
@ -151,7 +147,7 @@ class ElectronicStructure:
|
||||
|
||||
# Spin factor
|
||||
sp_fac = 2.0 if ns == 1 and self.nc_flag == False else 1.0
|
||||
|
||||
|
||||
if self.nc_flag == False:
|
||||
den_mat = np.zeros((ns, nproj, nproj), dtype=float)
|
||||
overlap = np.zeros((ns, nproj, nproj), dtype=float)
|
||||
@ -184,9 +180,9 @@ class ElectronicStructure:
|
||||
out += " "
|
||||
out += ''.join(map("{0:12.7f}".format, dov))
|
||||
print(out)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
else:
|
||||
print("!! WARNING !! Non Collinear Routine")
|
||||
den_mat = np.zeros((ns, nproj, nproj), dtype=float)
|
||||
|
@ -289,6 +289,8 @@ def ctrl_output(conf_pars, el_struct, ng):
|
||||
|
||||
* *nk*: number of `k`-points
|
||||
|
||||
* *nkibz*: number of `k`-points in IBZ
|
||||
|
||||
* *ns*: number of spin channels
|
||||
|
||||
* *nc_flag*: collinear/noncollinear case (False/True)
|
||||
@ -309,6 +311,7 @@ def ctrl_output(conf_pars, el_struct, ng):
|
||||
# Construct the header dictionary
|
||||
head_dict['ngroups'] = ng
|
||||
head_dict['nk'] = el_struct.kmesh['nktot']
|
||||
head_dict['nkibz'] = el_struct.kmesh['nkibz']
|
||||
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])
|
||||
|
@ -72,22 +72,26 @@ def is_vasp_running(vasp_pid):
|
||||
pid_exists = mpi.bcast(pid_exists)
|
||||
return pid_exists
|
||||
|
||||
|
||||
def get_dft_energy():
|
||||
"""
|
||||
Reads energy from the last line of OSZICAR.
|
||||
Reads DFT energy from the last line of Vasp's OSZICAR or from vasptriqs.h5
|
||||
"""
|
||||
with open('OSZICAR', 'r') as f:
|
||||
nextline = f.readline()
|
||||
while nextline.strip():
|
||||
line = nextline
|
||||
nextline = f.readline()
|
||||
# print "OSZICAR: ", line[:-1]
|
||||
h5_energy = False
|
||||
if os.path.isfile('vaspout.h5'):
|
||||
with HDFArchive('vaspout.h5', 'r') as h5:
|
||||
if 'oszicar' in h5['intermediate/ion_dynamics']:
|
||||
dft_energy = h5['intermediate/ion_dynamics/oszicar'][-1,1]
|
||||
h5_energy = True
|
||||
|
||||
try:
|
||||
# as backup use OSZICAR file
|
||||
if not h5_energy:
|
||||
with open('OSZICAR', 'r') as file:
|
||||
nextline = file.readline()
|
||||
while nextline.strip():
|
||||
line = nextline
|
||||
nextline = file.readline()
|
||||
dft_energy = float(line.split()[2])
|
||||
except ValueError:
|
||||
print("Cannot read energy from OSZICAR, setting it to zero")
|
||||
dft_energy = 0.0
|
||||
|
||||
return dft_energy
|
||||
|
||||
@ -98,9 +102,8 @@ class bcolors:
|
||||
YELLOW = '\033[93m'
|
||||
RED = '\033[91m'
|
||||
ENDC = '\033[0m'
|
||||
#
|
||||
|
||||
# Main self-consistent cycle
|
||||
#
|
||||
def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
|
||||
"""
|
||||
"""
|
||||
@ -117,15 +120,14 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
|
||||
mpi.barrier()
|
||||
while is_vasp_lock_present():
|
||||
time.sleep(1)
|
||||
# if debug: print bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC
|
||||
if debug: print(bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC)
|
||||
if not is_vasp_running(vasp_pid):
|
||||
mpi.report(" VASP stopped")
|
||||
vasp_running = False
|
||||
break
|
||||
|
||||
# Tell VASP to stop if the maximum number of iterations is reached
|
||||
|
||||
|
||||
# Tell VASP to stop if the maximum number of iterations is reached
|
||||
|
||||
if debug: print(bcolors.MAGENTA + "rank %s"%(mpi.rank) + bcolors.ENDC)
|
||||
err = 0
|
||||
exc = None
|
||||
@ -161,7 +163,7 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
|
||||
# electron.F around line 644
|
||||
iter_dft = 0
|
||||
|
||||
if vasp_version == 'standard':
|
||||
if vasp_version == 'standard' or vasp_version == 'ncl':
|
||||
copyfile(src='GAMMA',dst='GAMMA_recent')
|
||||
while iter_dft < n_iter_dft:
|
||||
# insert recalculation of GAMMA here
|
||||
@ -190,7 +192,7 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
|
||||
vasp_running = False
|
||||
break
|
||||
iter_dft += 1
|
||||
if vasp_version == 'standard':
|
||||
if vasp_version == 'standard' or vasp_version == 'ncl':
|
||||
copyfile(src='GAMMA_recent',dst='GAMMA')
|
||||
iter += 1
|
||||
if iter == n_iter:
|
||||
@ -253,8 +255,8 @@ def main():
|
||||
except KeyError:
|
||||
vasp_version = 'standard'
|
||||
|
||||
if vasp_version != 'standard' and vasp_version != 'no_gamma_write':
|
||||
raise Exception('vasp_version has to be standard or no_gamma_write')
|
||||
#if vasp_version != 'standard' and vasp_version != 'no_gamma_write':
|
||||
# raise Exception('vasp_version has to be standard or no_gamma_write')
|
||||
|
||||
# if len(sys.argv) > 1:
|
||||
# vasp_path = sys.argv[1]
|
||||
|
@ -1,4 +1,3 @@
|
||||
|
||||
################################################################################
|
||||
#
|
||||
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
||||
@ -40,9 +39,13 @@ r"""
|
||||
import logging
|
||||
import numpy as np
|
||||
import re
|
||||
import os
|
||||
import time
|
||||
from h5 import HDFArchive
|
||||
|
||||
log = logging.getLogger('plovasp.vaspio')
|
||||
|
||||
|
||||
def read_lines(filename):
|
||||
r"""
|
||||
Generator of lines for a file
|
||||
@ -56,6 +59,7 @@ def read_lines(filename):
|
||||
for line in f:
|
||||
yield line
|
||||
|
||||
|
||||
################################################################################
|
||||
################################################################################
|
||||
#
|
||||
@ -67,48 +71,63 @@ class VaspData:
|
||||
"""
|
||||
Container class for all VASP data.
|
||||
"""
|
||||
|
||||
def __init__(self, vasp_dir, read_all=True, efermi_required=True):
|
||||
self.vasp_dir = vasp_dir
|
||||
|
||||
self.plocar = Plocar()
|
||||
self.poscar = Poscar()
|
||||
self.kpoints = Kpoints()
|
||||
self.eigenval = Eigenval()
|
||||
self.doscar = Doscar()
|
||||
# read from vaspout.h5 if possible
|
||||
vasph5 = os.path.isfile(os.path.join(vasp_dir, 'vaspout.h5'))
|
||||
if vasph5:
|
||||
log.warning("Reading from vaspout.h5")
|
||||
h5path = os.path.join(vasp_dir, 'vaspout.h5')
|
||||
# give VASP some time to write the file
|
||||
with HDFArchive(h5path, 'r') as archive:
|
||||
if 'locproj' not in archive['results']:
|
||||
time.sleep(2)
|
||||
self.plocar = h5Plocar(h5path)
|
||||
self.poscar = h5Poscar(h5path)
|
||||
self.kpoints = h5Kpoints(h5path)
|
||||
self.eigenval = h5Eigenval(h5path, self.kpoints.ksymmap)
|
||||
self.doscar = h5Doscar(h5path)
|
||||
else:
|
||||
self.plocar = Plocar()
|
||||
self.poscar = Poscar()
|
||||
self.kpoints = Kpoints()
|
||||
self.eigenval = Eigenval()
|
||||
self.doscar = Doscar()
|
||||
|
||||
if read_all:
|
||||
self.plocar.from_file(vasp_dir)
|
||||
self.poscar.from_file(vasp_dir)
|
||||
self.kpoints.from_file(vasp_dir)
|
||||
try:
|
||||
self.eigenval.from_file(vasp_dir)
|
||||
except (IOError, StopIteration):
|
||||
self.eigenval.eigs = None
|
||||
self.eigenval.ferw = None
|
||||
log.warning("Error reading from EIGENVAL, trying LOCPROJ...")
|
||||
if read_all:
|
||||
self.plocar.from_file(vasp_dir)
|
||||
self.poscar.from_file(vasp_dir)
|
||||
self.kpoints.from_file(vasp_dir)
|
||||
try:
|
||||
self.eigenval.from_file(vasp_dir)
|
||||
except (IOError, StopIteration):
|
||||
self.eigenval.eigs = None
|
||||
self.eigenval.ferw = None
|
||||
log.warning("Error reading from EIGENVAL, trying LOCPROJ...")
|
||||
|
||||
try:
|
||||
self.doscar.from_file(vasp_dir)
|
||||
except (IOError, StopIteration):
|
||||
if efermi_required:
|
||||
log.warning("Error reading Efermi from DOSCAR, trying LOCPROJ...")
|
||||
try:
|
||||
self.plocar.efermi
|
||||
self.doscar.efermi = self.plocar.efermi
|
||||
except NameError:
|
||||
raise Exception("Efermi cannot be read from DOSCAR or LOCPROJ")
|
||||
else:
|
||||
# TODO: This a hack. Find out a way to determine ncdij without DOSCAR
|
||||
log.warning("Error reading Efermi from DOSCAR, taking from config")
|
||||
self.doscar.ncdij = self.plocar.nspin
|
||||
try:
|
||||
self.doscar.from_file(vasp_dir)
|
||||
except (IOError, StopIteration):
|
||||
if efermi_required:
|
||||
log.warning("Error reading Efermi from DOSCAR, trying LOCPROJ...")
|
||||
try:
|
||||
self.plocar.efermi
|
||||
self.doscar.efermi = self.plocar.efermi
|
||||
except NameError:
|
||||
raise Exception("Efermi cannot be read from DOSCAR or LOCPROJ")
|
||||
else:
|
||||
# TODO: This a hack. Find out a way to determine ncdij without DOSCAR
|
||||
log.warning("Error reading Efermi from DOSCAR, taking from config")
|
||||
self.doscar.ncdij = self.plocar.nspin
|
||||
|
||||
################################################################################
|
||||
################################################################################
|
||||
|
||||
##########################
|
||||
#
|
||||
# class Plocar
|
||||
#
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
class Plocar:
|
||||
"""
|
||||
Class containing raw PLO data from VASP.
|
||||
@ -119,6 +138,7 @@ class Plocar:
|
||||
- *ferw* (array(nion, ns, nk, nb)) : Fermi weights from VASP
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
self.plo = None
|
||||
self.proj_params = None
|
||||
@ -135,15 +155,14 @@ class Plocar:
|
||||
plocar_filename (str) : filename [default = 'PLOCAR']
|
||||
|
||||
"""
|
||||
# Add a slash to the path name if necessary
|
||||
# Add a slash to the path name if necessary
|
||||
if vasp_dir[-1] != '/':
|
||||
vasp_dir += '/'
|
||||
|
||||
# self.params, self.plo, self.ferw = c_plocar_io.read_plocar(vasp_dir + plocar_filename)
|
||||
# self.proj_params, self.plo = self.temp_parser(projcar_filename=vasp_dir + "PROJCAR", locproj_filename=vasp_dir + "LOCPROJ")
|
||||
# self.params, self.plo, self.ferw = c_plocar_io.read_plocar(vasp_dir + plocar_filename)
|
||||
# self.proj_params, self.plo = self.temp_parser(projcar_filename=vasp_dir + "PROJCAR", locproj_filename=vasp_dir + "LOCPROJ")
|
||||
self.proj_params, self.plo = self.locproj_parser(locproj_filename=vasp_dir + "LOCPROJ")
|
||||
|
||||
|
||||
def locproj_parser(self, locproj_filename='LOCPROJ'):
|
||||
r"""
|
||||
Parses LOCPROJ (for VASP >= 5.4.2) to get VASP projectors.
|
||||
@ -159,20 +178,20 @@ class Plocar:
|
||||
|
||||
def lm_to_l_m(lm):
|
||||
l = int(np.sqrt(lm))
|
||||
m = lm - l*l
|
||||
m = lm - l * l
|
||||
return l, m
|
||||
|
||||
# Read the first line of LOCPROJ to get the dimensions
|
||||
# Read the first line of LOCPROJ to get the dimensions
|
||||
with open(locproj_filename, 'rt') as f:
|
||||
line = f.readline()
|
||||
line = line.split("#")[0]
|
||||
sline = line.split()
|
||||
self.ncdij, nk, self.nband, nproj = list(map(int, sline[0:4]))
|
||||
|
||||
|
||||
# VASP.6.
|
||||
self.nspin = self.ncdij if self.ncdij < 4 else 1
|
||||
log.debug("ISPIN is {}".format(self.nspin))
|
||||
|
||||
|
||||
self.nspin_band = 2 if self.ncdij == 2 else 1
|
||||
|
||||
try:
|
||||
@ -185,7 +204,7 @@ class Plocar:
|
||||
|
||||
iproj_site = 0
|
||||
is_first_read = True
|
||||
|
||||
|
||||
# VASP.6.
|
||||
if self.ncdij == 4:
|
||||
self.nc_flag = 1
|
||||
@ -195,7 +214,7 @@ class Plocar:
|
||||
|
||||
log.debug("NC FLAG : {}".format(self.nc_flag))
|
||||
|
||||
# First read the header block with orbital labels
|
||||
# First read the header block with orbital labels
|
||||
line = self.search_for(f, "^ *ISITE")
|
||||
ip = 0
|
||||
while line:
|
||||
@ -204,23 +223,23 @@ class Plocar:
|
||||
label = sline[-1].strip()
|
||||
lm = orb_labels.index(label)
|
||||
l, m = lm_to_l_m(lm)
|
||||
# ip_new = iproj_site * norb + il
|
||||
# ip_prev = (iproj_site - 1) * norb + il
|
||||
# ip_new = iproj_site * norb + il
|
||||
# ip_prev = (iproj_site - 1) * norb + il
|
||||
proj_params[ip]['label'] = label
|
||||
proj_params[ip]['isite'] = isite
|
||||
proj_params[ip]['l'] = l
|
||||
if self.nc_flag == True:
|
||||
if (ip % 2) == 0:
|
||||
proj_params[ip]['m'] = 2*m
|
||||
proj_params[ip]['m'] = 2 * m
|
||||
else:
|
||||
proj_params[ip]['m'] = 2*m + 1
|
||||
proj_params[ip]['m'] = 2 * m + 1
|
||||
else:
|
||||
proj_params[ip]['m'] = m
|
||||
|
||||
ip +=1
|
||||
|
||||
ip += 1
|
||||
|
||||
line = f.readline().strip()
|
||||
|
||||
|
||||
assert ip == nproj, "Number of projectors in the header is wrong in LOCPROJ"
|
||||
|
||||
self.eigs = np.zeros((nk, self.nband, self.nspin_band))
|
||||
@ -243,7 +262,7 @@ class Plocar:
|
||||
line = f.readline()
|
||||
sline = line.split()
|
||||
ctmp = complex(float(sline[1]), float(sline[2]))
|
||||
plo[ip, ispin, ik, ib] = ctmp
|
||||
plo[ip, ispin, ik, ib] = ctmp
|
||||
|
||||
print("Read parameters: LOCPROJ")
|
||||
for il, par in enumerate(proj_params):
|
||||
@ -251,7 +270,6 @@ class Plocar:
|
||||
|
||||
return proj_params, plo
|
||||
|
||||
|
||||
def search_for(self, f, patt):
|
||||
r"""
|
||||
Reads file 'f' until pattern 'patt' is encountered and returns
|
||||
@ -265,13 +283,11 @@ class Plocar:
|
||||
return line
|
||||
|
||||
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
#
|
||||
# class Poscar
|
||||
#
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
class Poscar:
|
||||
"""
|
||||
Class containing POSCAR data from VASP.
|
||||
@ -285,6 +301,7 @@ class Poscar:
|
||||
- q_types ([numpy.array((nions, 3), dtype=float)]) : a list of
|
||||
arrays each containing fractional coordinates of ions of a given type
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
self.q_cart = None
|
||||
|
||||
@ -299,68 +316,69 @@ class Poscar:
|
||||
plocar_filename (str) : filename [default = 'POSCAR']
|
||||
|
||||
"""
|
||||
# Convenince local function
|
||||
|
||||
# Convenince local function
|
||||
def readline_remove_comments():
|
||||
return next(f).split('!')[0].split('#')[0].strip()
|
||||
|
||||
# Add a slash to the path name if necessary
|
||||
# Add a slash to the path name if necessary
|
||||
if vasp_dir[-1] != '/':
|
||||
vasp_dir += '/'
|
||||
|
||||
f = read_lines(vasp_dir + poscar_filename)
|
||||
# Comment line
|
||||
# Comment line
|
||||
comment = next(f).rstrip()
|
||||
print(" Found POSCAR, title line: %s"%(comment))
|
||||
print(" Found POSCAR, title line: %s" % (comment))
|
||||
|
||||
# Read scale
|
||||
# Read scale
|
||||
sline = readline_remove_comments()
|
||||
ascale = float(sline)
|
||||
# Read lattice vectors
|
||||
# Read lattice vectors
|
||||
self.a_brav = np.zeros((3, 3))
|
||||
for ia in range(3):
|
||||
sline = readline_remove_comments()
|
||||
self.a_brav[ia, :] = list(map(float, sline.split()))
|
||||
# Negative scale means that it is a volume scale
|
||||
# Negative scale means that it is a volume scale
|
||||
if ascale < 0:
|
||||
vscale = -ascale
|
||||
vol = np.linalg.det(self.a_brav)
|
||||
ascale = (vscale / vol)**(1.0/3)
|
||||
ascale = (vscale / vol) ** (1.0 / 3)
|
||||
|
||||
self.a_brav *= ascale
|
||||
|
||||
# Depending on the version of VASP there could be
|
||||
# an extra line with element names
|
||||
# Depending on the version of VASP there could be
|
||||
# an extra line with element names
|
||||
sline = readline_remove_comments()
|
||||
try:
|
||||
# Old v4.6 format: no element names
|
||||
# Old v4.6 format: no element names
|
||||
self.nions = list(map(int, sline.split()))
|
||||
self.el_names = ['El%i'%(i) for i in range(len(self.nions))]
|
||||
self.el_names = ['El%i' % (i) for i in range(len(self.nions))]
|
||||
except ValueError:
|
||||
# New v5.x format: read element names first
|
||||
# New v5.x format: read element names first
|
||||
self.el_names = sline.split()
|
||||
sline = readline_remove_comments()
|
||||
self.nions = list(map(int, sline.split()))
|
||||
|
||||
# Set the number of atom sorts (types) and the total
|
||||
# number of atoms in the unit cell
|
||||
# Set the number of atom sorts (types) and the total
|
||||
# number of atoms in the unit cell
|
||||
self.ntypes = len(self.nions)
|
||||
self.nq = sum(self.nions)
|
||||
|
||||
# Check for the line 'Selective dynamics' (and ignore it)
|
||||
# Check for the line 'Selective dynamics' (and ignore it)
|
||||
sline = readline_remove_comments()
|
||||
if sline[0].lower() == 's':
|
||||
sline = readline_remove_comments()
|
||||
|
||||
# Check whether coordinates are cartesian or fractional
|
||||
# Check whether coordinates are cartesian or fractional
|
||||
cartesian = (sline[0].lower() in 'ck')
|
||||
# determine reciprocal basis in units of 2*pi
|
||||
self.kpt_basis = np.linalg.inv(self.a_brav.T)
|
||||
|
||||
# Read atomic positions
|
||||
# Read atomic positions
|
||||
self.q_types = []
|
||||
self.type_of_ion = []
|
||||
for it in range(self.ntypes):
|
||||
# Array mapping ion index to type
|
||||
# Array mapping ion index to type
|
||||
self.type_of_ion += self.nions[it] * [it]
|
||||
|
||||
q_at_it = np.zeros((self.nions[it], 3))
|
||||
@ -377,38 +395,35 @@ class Poscar:
|
||||
print(" Number of types:", self.ntypes)
|
||||
print(" Number of ions for each type:", self.nions)
|
||||
|
||||
# print
|
||||
# print " Coords:"
|
||||
# for it in range(ntypes):
|
||||
# print " Element:", el_names[it]
|
||||
# print q_at[it]
|
||||
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
#
|
||||
# class Kpoints
|
||||
#
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
class Kpoints:
|
||||
"""
|
||||
Class describing k-points and optionally tetrahedra.
|
||||
|
||||
Properties:
|
||||
- nktot (int) : total number of k-points in the IBZ
|
||||
- nktot (int) : total number of k-points in the BZ
|
||||
- nkibz (int) : number of k-points in the IBZ
|
||||
- kpts (numpy.array((nktot, 3), dtype=float)) : k-point vectors (fractional coordinates)
|
||||
- ntet (int) : total number of k-point tetrahedra
|
||||
- itet (numpy.array((ntet, 5), dtype=float) : array of tetrahedra
|
||||
- volt (float) : volume of a tetrahedron (the k-grid is assumed to
|
||||
be uniform)
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
self.kpts = None
|
||||
self.nktot = None
|
||||
self.nkibz = None
|
||||
self.kwghts = None
|
||||
#
|
||||
# Reads IBZKPT file
|
||||
#
|
||||
|
||||
#
|
||||
# Reads IBZKPT file
|
||||
#
|
||||
def from_file(self, vasp_dir='./', ibz_filename='IBZKPT'):
|
||||
r"""
|
||||
Reads from IBZKPT: k-points and optionally
|
||||
@ -422,17 +437,20 @@ class Kpoints:
|
||||
|
||||
"""
|
||||
|
||||
# Add a slash to the path name if necessary
|
||||
# Add a slash to the path name if necessary
|
||||
if vasp_dir[-1] != '/':
|
||||
vasp_dir += '/'
|
||||
|
||||
ibz_file = read_lines(vasp_dir + ibz_filename)
|
||||
|
||||
# Skip comment line
|
||||
# Skip comment line
|
||||
line = next(ibz_file)
|
||||
# Number of k-points
|
||||
# Number of k-points
|
||||
line = next(ibz_file)
|
||||
self.nktot = int(line.strip().split()[0])
|
||||
# when reading from IBZKPT file we do not know the full number of k-points, i.e.
|
||||
# works only ISYM=-1
|
||||
self.nkibz = self.nktot
|
||||
|
||||
print()
|
||||
print(" {0:>26} {1:d}".format("Total number of k-points:", self.nktot))
|
||||
@ -440,7 +458,7 @@ class Kpoints:
|
||||
self.kpts = np.zeros((self.nktot, 3))
|
||||
self.kwghts = np.zeros((self.nktot))
|
||||
|
||||
# Skip comment line
|
||||
# Skip comment line
|
||||
line = next(ibz_file)
|
||||
for ik in range(self.nktot):
|
||||
line = next(ibz_file)
|
||||
@ -450,12 +468,12 @@ class Kpoints:
|
||||
|
||||
self.kwghts /= self.nktot
|
||||
|
||||
# Attempt to read tetrahedra
|
||||
# Skip comment line ("Tetrahedra")
|
||||
# Attempt to read tetrahedra
|
||||
# Skip comment line ("Tetrahedra")
|
||||
try:
|
||||
line = next(ibz_file)
|
||||
|
||||
# Number of tetrahedra and volume = 1/(6*nkx*nky*nkz)
|
||||
# Number of tetrahedra and volume = 1/(6*nkx*nky*nkz)
|
||||
line = next(ibz_file)
|
||||
sline = line.split()
|
||||
self.ntet = int(sline[0])
|
||||
@ -463,35 +481,26 @@ class Kpoints:
|
||||
|
||||
print(" {0:>26} {1:d}".format("Total number of tetrahedra:", self.ntet))
|
||||
|
||||
# Traditionally, itet[it, 0] contains multiplicity
|
||||
# Traditionally, itet[it, 0] contains multiplicity
|
||||
self.itet = np.zeros((self.ntet, 5), dtype=int)
|
||||
for it in range(self.ntet):
|
||||
line = next(ibz_file)
|
||||
self.itet[it, :] = list(map(int, line.split()[:5]))
|
||||
line = next(ibz_file)
|
||||
self.itet[it, :] = list(map(int, line.split()[:5]))
|
||||
except StopIteration as ValueError:
|
||||
print(" No tetrahedron data found in %s. Skipping..."%(ibz_filename))
|
||||
print(" No tetrahedron data found in %s. Skipping..." % (ibz_filename))
|
||||
self.ntet = 0
|
||||
|
||||
# data = { 'nktot': nktot,
|
||||
# 'kpts': kpts,
|
||||
# 'ntet': ntet,
|
||||
# 'itet': itet,
|
||||
# 'volt': volt }
|
||||
#
|
||||
# return data
|
||||
|
||||
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
#
|
||||
# class Eigenval
|
||||
#
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
class Eigenval:
|
||||
"""
|
||||
Class containing Kohn-Sham-eigenvalues data from VASP (EIGENVAL file).
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
self.eigs = None
|
||||
self.ferw = None
|
||||
@ -503,43 +512,43 @@ class Eigenval:
|
||||
then used to check the consistency of files read.
|
||||
"""
|
||||
|
||||
# Add a slash to the path name if necessary
|
||||
# Add a slash to the path name if necessary
|
||||
if vasp_dir[-1] != '/':
|
||||
vasp_dir += '/'
|
||||
|
||||
f = read_lines(vasp_dir + eig_filename)
|
||||
|
||||
# First line: only the first and the last number out of four
|
||||
# are used; these are 'nions' and 'ispin'
|
||||
# First line: only the first and the last number out of four
|
||||
# are used; these are 'nions' and 'ispin'
|
||||
sline = next(f).split()
|
||||
self.nq = int(sline[0])
|
||||
self.ispin = int(sline[3])
|
||||
|
||||
# Second line: cell volume and lengths of lattice vectors (skip)
|
||||
# Second line: cell volume and lengths of lattice vectors (skip)
|
||||
sline = next(f)
|
||||
|
||||
# Third line: temperature (skip)
|
||||
# Third line: temperature (skip)
|
||||
sline = next(f)
|
||||
|
||||
# Fourth and fifth line: useless
|
||||
# Fourth and fifth line: useless
|
||||
sline = next(f)
|
||||
sline = next(f)
|
||||
|
||||
# Sixth line: NELECT, NKTOT, NBTOT
|
||||
# Sixth line: NELECT, NKTOT, NBTOT
|
||||
sline = next(f).split()
|
||||
self.nelect = int(sline[0])
|
||||
self.nktot = int(sline[1])
|
||||
self.nband = int(sline[2])
|
||||
|
||||
# Set of eigenvalues and k-points
|
||||
# Set of eigenvalues and k-points
|
||||
self.kpts = np.zeros((self.nktot, 3))
|
||||
self.kwghts = np.zeros((self.nktot,))
|
||||
self.eigs = np.zeros((self.nktot, self.nband, self.ispin))
|
||||
self.ferw = np.zeros((self.nktot, self.nband, self.ispin))
|
||||
|
||||
for ik in range(self.nktot):
|
||||
sline = next(f) # Empty line
|
||||
sline = next(f) # k-point info
|
||||
sline = next(f) # Empty line
|
||||
sline = next(f) # k-point info
|
||||
tmp = list(map(float, sline.split()))
|
||||
self.kpts[ik, :] = tmp[:3]
|
||||
self.kwghts[ik] = tmp[3]
|
||||
@ -548,21 +557,20 @@ class Eigenval:
|
||||
sline = next(f).split()
|
||||
tmp = list(map(float, sline))
|
||||
assert len(tmp) == 2 * self.ispin + 1, "EIGENVAL file is incorrect (probably from old versions of VASP)"
|
||||
self.eigs[ik, ib, :] = tmp[1:self.ispin+1]
|
||||
self.ferw[ik, ib, :] = tmp[self.ispin+1:]
|
||||
self.eigs[ik, ib, :] = tmp[1:self.ispin + 1]
|
||||
self.ferw[ik, ib, :] = tmp[self.ispin + 1:]
|
||||
|
||||
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
#
|
||||
# class Doscar
|
||||
#
|
||||
################################################################################
|
||||
################################################################################
|
||||
##########################
|
||||
class Doscar:
|
||||
"""
|
||||
Class containing some data from DOSCAR
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
self.ncdij = None
|
||||
self.efermi = None
|
||||
@ -572,38 +580,40 @@ class Doscar:
|
||||
Reads only E_Fermi from DOSCAR.
|
||||
"""
|
||||
|
||||
# Add a slash to the path name if necessary
|
||||
# Add a slash to the path name if necessary
|
||||
if vasp_dir[-1] != '/':
|
||||
vasp_dir += '/'
|
||||
|
||||
f = read_lines(vasp_dir + dos_filename)
|
||||
|
||||
# First line: NION, NION, JOBPAR, NCDIJ
|
||||
# First line: NION, NION, JOBPAR, NCDIJ
|
||||
sline = next(f).split()
|
||||
|
||||
# Skip next 4 lines
|
||||
|
||||
# Skip next 4 lines
|
||||
for _ in range(4):
|
||||
sline = next(f)
|
||||
|
||||
# Sixth line: EMAX, EMIN, NEDOS, EFERMI, 1.0
|
||||
# Sixth line: EMAX, EMIN, NEDOS, EFERMI, 1.0
|
||||
sline = next(f).split()
|
||||
self.efermi = float(sline[3])
|
||||
|
||||
|
||||
# TODO: implement output of SYMMCAR in VASP and read it here
|
||||
################################################################
|
||||
##########################
|
||||
#
|
||||
# Reads SYMMCAR
|
||||
#
|
||||
################################################################
|
||||
##########################
|
||||
def read_symmcar(vasp_dir, symm_filename='SYMMCAR'):
|
||||
"""
|
||||
Reads SYMMCAR.
|
||||
"""
|
||||
# Shorthand for simple parsing
|
||||
|
||||
# Shorthand for simple parsing
|
||||
def extract_int_par(parname):
|
||||
return int(re.findall(parname + '\s*=\s*(\d+)', line)[-1])
|
||||
|
||||
# Add a slash to the path name if necessary
|
||||
# Add a slash to the path name if necessary
|
||||
if vasp_dir[-1] != '/':
|
||||
vasp_dir += '/'
|
||||
|
||||
@ -614,11 +624,11 @@ def read_symmcar(vasp_dir, symm_filename='SYMMCAR'):
|
||||
|
||||
line = next(sym_file)
|
||||
ntrans = extract_int_par('NPCELL')
|
||||
# Lmax
|
||||
# Lmax
|
||||
line = next(sym_file)
|
||||
lmax = extract_int_par('LMAX')
|
||||
mmax = 2 * lmax + 1
|
||||
# Nion
|
||||
# Nion
|
||||
line = next(sym_file)
|
||||
nion = extract_int_par('NION')
|
||||
|
||||
@ -627,22 +637,22 @@ def read_symmcar(vasp_dir, symm_filename='SYMMCAR'):
|
||||
print(" {0:>26} {1:d}".format("Number of ions:", nion))
|
||||
print(" {0:>26} {1:d}".format("L_max:", lmax))
|
||||
|
||||
rot_mats = np.zeros((nrot, lmax+1, mmax, mmax))
|
||||
rot_mats = np.zeros((nrot, lmax + 1, mmax, mmax))
|
||||
rot_map = np.zeros((nrot, ntrans, nion), dtype=int)
|
||||
|
||||
for irot in range(nrot):
|
||||
# Empty line
|
||||
# Empty line
|
||||
line = next(sym_file)
|
||||
# IROT index (skip it)
|
||||
# IROT index (skip it)
|
||||
line = next(sym_file)
|
||||
# ISYMOP matrix (can be also skipped)
|
||||
# ISYMOP matrix (can be also skipped)
|
||||
line = next(sym_file)
|
||||
line = next(sym_file)
|
||||
line = next(sym_file)
|
||||
|
||||
# Skip comment " Permutation map..."
|
||||
# Skip comment " Permutation map..."
|
||||
line = next(sym_file)
|
||||
# Permutations (in chunks of 20 indices per line)
|
||||
# Permutations (in chunks of 20 indices per line)
|
||||
for it in range(ntrans):
|
||||
for ibl in range((nion - 1) // 20 + 1):
|
||||
i1 = ibl * 20
|
||||
@ -652,12 +662,169 @@ def read_symmcar(vasp_dir, symm_filename='SYMMCAR'):
|
||||
|
||||
for l in range(lmax + 1):
|
||||
mmax = 2 * l + 1
|
||||
# Comment: "L = ..."
|
||||
# Comment: "L = ..."
|
||||
line = next(sym_file)
|
||||
for m in range(mmax):
|
||||
line = next(sym_file)
|
||||
rot_mats[irot, l, m, :mmax] = list(map(float, line.split()[:mmax]))
|
||||
|
||||
data.update({ 'nrot': nrot, 'ntrans': ntrans,
|
||||
'lmax': lmax, 'nion': nion,
|
||||
'sym_rots': rot_mats, 'perm_map': rot_map })
|
||||
data.update({'nrot': nrot, 'ntrans': ntrans,
|
||||
'lmax': lmax, 'nion': nion,
|
||||
'sym_rots': rot_mats, 'perm_map': rot_map})
|
||||
|
||||
|
||||
class h5Poscar:
|
||||
|
||||
def __init__(self, h5path):
|
||||
# self.q_cart = None
|
||||
|
||||
with HDFArchive(h5path, 'a') as archive:
|
||||
struct = archive['results/positions']
|
||||
ascale = struct['scale']
|
||||
self.a_brav = struct['lattice_vectors']
|
||||
self.nions = struct['number_ion_types']
|
||||
self.el_names = struct['ion_types']
|
||||
direct = struct['direct_coordinates']
|
||||
qcoord = struct['position_ions']
|
||||
|
||||
if ascale < 0:
|
||||
vscale = -ascale
|
||||
vol = np.linalg.det(self.a_brav)
|
||||
ascale = (vscale / vol) ** (1.0 / 3)
|
||||
self.a_brav *= ascale
|
||||
# determine reciprocal basis in units of 2*pi
|
||||
self.kpt_basis = np.linalg.inv(self.a_brav.T)
|
||||
|
||||
# Set the number of atom sorts (types) and the total
|
||||
# number of atoms in the unit cell
|
||||
self.ntypes = len(self.nions)
|
||||
self.nq = sum(self.nions)
|
||||
|
||||
# Read atomic positions
|
||||
self.q_types = []
|
||||
self.type_of_ion = []
|
||||
for it in range(self.ntypes):
|
||||
# Array mapping ion index to type
|
||||
self.type_of_ion += self.nions[it] * [it]
|
||||
q_at_it = np.zeros((self.nions[it], 3))
|
||||
for iq in range(self.nions[it]):
|
||||
if not direct:
|
||||
qcoord = np.dot(self.kpt_basis, qcoord[iq + it])
|
||||
q_at_it[iq, :] = qcoord[iq + it]
|
||||
|
||||
self.q_types.append(q_at_it)
|
||||
|
||||
print(" Total number of ions:", self.nq)
|
||||
print(" Number of types:", self.ntypes)
|
||||
print(" Number of ions for each type:", self.nions)
|
||||
|
||||
|
||||
class h5Kpoints:
|
||||
|
||||
def __init__(self, h5path):
|
||||
|
||||
# h5path = './vasptriqs.h5'
|
||||
with HDFArchive(h5path, 'a') as archive:
|
||||
kpoints = archive['results/electron_eigenvalues']
|
||||
self.nkibz = kpoints['kpoints']
|
||||
self.kpts = kpoints['kpoint_coords_full']
|
||||
self.nktot = len(self.kpts)
|
||||
self.kwghts = kpoints['kpoints_symmetry_weight_full']
|
||||
self.ksymmap = kpoints['kpoints_symmetry_mapping']
|
||||
self.ksymmap -= 1
|
||||
try:
|
||||
self.ntet = kpoints['num_tetrahedra']
|
||||
self.vtet = kpoints['volume_weight_tetrahedra']
|
||||
self.itet = kpoints['coordinate_id_tetrahedra']
|
||||
except KeyError:
|
||||
print(" No tetrahedron data found in vaspout.h5. Skipping...")
|
||||
self.ntet = 0
|
||||
|
||||
print()
|
||||
print(" {0:>26} {1:d}".format("Reduced number of k-points:", self.nkibz))
|
||||
print(" {0:>26} {1:d}".format("Total number of k-points:", self.nktot))
|
||||
print(" {0:>26} {1:d}".format("Total number of tetrahedra:", self.ntet))
|
||||
|
||||
|
||||
class h5Eigenval:
|
||||
|
||||
def __init__(self, h5path, symmap):
|
||||
with HDFArchive(h5path, 'a') as archive:
|
||||
self.eigs = archive['results/electron_eigenvalues']['eigenvalues']
|
||||
self.eigs = self.eigs[:, symmap, :]
|
||||
self.ferw = archive['results/electron_eigenvalues']['fermiweights']
|
||||
self.ferw = self.ferw[:, symmap, :]
|
||||
# TODO Change the format in VASP to have [kpoints, bands, spin]
|
||||
self.eigs = np.transpose(self.eigs, (1, 2, 0))
|
||||
self.ferw = np.transpose(self.ferw, (1, 2, 0))
|
||||
|
||||
|
||||
class h5Doscar:
|
||||
|
||||
def __init__(self, h5path):
|
||||
with HDFArchive(h5path, 'a') as archive:
|
||||
self.efermi = archive['results/electron_dos']['efermi']
|
||||
|
||||
|
||||
class h5Plocar():
|
||||
|
||||
def __init__(self, h5path):
|
||||
with HDFArchive(h5path, 'a') as archive:
|
||||
plo = np.array(archive['results/locproj']['data'])
|
||||
self.nc_flag = int(archive['results/locproj/parameters']['lnoncollinear'])
|
||||
|
||||
self.nproj = plo.shape[0]
|
||||
self.ncdij = plo.shape[1]
|
||||
nk = plo.shape[2]
|
||||
self.nband = plo.shape[3]
|
||||
|
||||
self.nspin = self.ncdij if self.ncdij < 4 else 1
|
||||
self.nspin_band = 2 if self.ncdij == 2 else 1
|
||||
|
||||
log.debug("ISPIN is {}".format(self.nspin))
|
||||
log.debug("NC FLAG : {}".format(self.nc_flag))
|
||||
|
||||
# self.proj_params, self.plo = self.locproj_parser(locproj_filename=vasp_dir + "LOCPROJ")
|
||||
|
||||
orb_labels = ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2",
|
||||
"fy(3x2-y2)", "fxyz", "fyz2", "fz3", "fxz2", "fz(x2-y2)", "fx(x2-3y2)"]
|
||||
|
||||
self.plo = np.zeros((self.nproj, self.nspin, nk, self.nband), dtype=complex)
|
||||
|
||||
for proj in range(self.nproj):
|
||||
for spin in range(self.nspin):
|
||||
for kpt in range(nk):
|
||||
for band in range(self.nband):
|
||||
real_plo = plo[proj, spin, kpt, band, 0] # Real part
|
||||
imag_plo = plo[proj, spin, kpt, band, 1] # Imaginary part
|
||||
self.plo[proj, spin, kpt, band] = complex(real_plo, imag_plo)
|
||||
|
||||
def lm_to_l_m(lm):
|
||||
l = int(np.sqrt(lm))
|
||||
m = lm - l * l
|
||||
return l, m
|
||||
|
||||
self.proj_params = [{} for i in range(self.nproj)]
|
||||
with HDFArchive(h5path, 'a') as archive:
|
||||
for it in range(self.nproj):
|
||||
projectors = archive['results/locproj']['parameters']
|
||||
self.proj_params[it]['label'] = projectors['ang_type'][it]
|
||||
self.proj_params[it]['isite'] = projectors['site'][it]
|
||||
self.proj_params[it]['coord'] = projectors['coordinates'][it]
|
||||
|
||||
for it in range(self.nproj):
|
||||
lm = orb_labels.index(self.proj_params[it]['label'].strip())
|
||||
l, m = lm_to_l_m(lm)
|
||||
self.proj_params[it]['l'] = l
|
||||
if self.nc_flag == True:
|
||||
if (it % 2) == 0:
|
||||
self.proj_params[it]['m'] = 2 * m
|
||||
else:
|
||||
self.proj_params[it]['m'] = 2 * m + 1
|
||||
else:
|
||||
self.proj_params[it]['m'] = m
|
||||
# assert ip == nproj, "Number of projectors in the header is wrong in LOCPROJ"
|
||||
|
||||
print("Read parameters: LOCPROJ")
|
||||
for il, par in enumerate(self.proj_params):
|
||||
print(il, " -> ", par)
|
||||
|
@ -158,6 +158,8 @@ class VaspConverter(ConverterTools):
|
||||
|
||||
ng = ctrl_head['ngroups']
|
||||
n_k = ctrl_head['nk']
|
||||
n_k_ibz = ctrl_head['nkibz']
|
||||
|
||||
# Note the difference in name conventions!
|
||||
SP = ctrl_head['ns'] - 1
|
||||
SO = ctrl_head['nc_flag']
|
||||
@ -387,7 +389,7 @@ class VaspConverter(ConverterTools):
|
||||
with HDFArchive(self.hdf_file,'a') as ar:
|
||||
if not (self.dft_subgrp in ar): ar.create_group(self.dft_subgrp)
|
||||
# The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten!
|
||||
things_to_save = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required',
|
||||
things_to_save = ['energy_unit','n_k', 'k_dep_projection','SP','SO','charge_below','density_required',
|
||||
'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat',
|
||||
'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights',
|
||||
'hopping','n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr','proj_or_hk',
|
||||
@ -401,6 +403,8 @@ class VaspConverter(ConverterTools):
|
||||
ar[self.misc_subgrp]['dft_fermi_weights'] = f_weights
|
||||
ar[self.misc_subgrp]['kpts_cart'] = kpts_cart
|
||||
ar[self.misc_subgrp]['band_window'] = band_window
|
||||
if n_k_ibz is not None:
|
||||
ar[self.misc_subgrp]['n_k_ibz'] = n_k_ibz
|
||||
|
||||
# Symmetries are used, so now convert symmetry information for *correlated* orbitals:
|
||||
self.convert_symmetry_input(ctrl_head, orbits=self.corr_shells, symm_subgrp=self.symmcorr_subgrp)
|
||||
|
@ -477,8 +477,7 @@ class SumkDFT(object):
|
||||
gf_rotated.from_L_G_R(rot_mat[ish].conjugate(
|
||||
), gf_rotated, rot_mat[ish].transpose())
|
||||
else:
|
||||
gf_rotated.from_L_G_R(rot_mat[ish], gf_rotated, rot_mat[
|
||||
ish].conjugate().transpose())
|
||||
gf_rotated.from_L_G_R(rot_mat[ish], gf_rotated, rot_mat[ish].conjugate().transpose())
|
||||
|
||||
elif direction == 'toLocal':
|
||||
|
||||
@ -541,8 +540,7 @@ class SumkDFT(object):
|
||||
else: # Check that existing GF is consistent
|
||||
G_latt = self.G_latt
|
||||
GFsize = [gf.target_shape[0] for bname, gf in G_latt]
|
||||
unchangedsize = all([self.n_orbitals[ik, ntoi[spn[isp]]] == GFsize[
|
||||
isp] for isp in range(self.n_spin_blocks[self.SO])])
|
||||
unchangedsize = all([self.n_orbitals[ik, ntoi[spn[isp]]] == GFsize[isp] for isp in range(self.n_spin_blocks[self.SO])])
|
||||
if (not mesh is None) or (not unchangedsize):
|
||||
set_up_G_latt = True
|
||||
|
||||
@ -1564,7 +1562,7 @@ class SumkDFT(object):
|
||||
warn("WARNING: density_matrix: method 'using_point_integration' is deprecated. Use 'density_matrix_using_point_integration' instead. All additionally provided arguments are ignored.")
|
||||
dens_mat = self.density_matrix_using_point_integration()
|
||||
else:
|
||||
raise ValueError("density_matrix: the method '%s' is not supported." % method)
|
||||
raise ValueError("density_matrix: the method '%s' is not supported." % method)
|
||||
|
||||
return dens_mat
|
||||
|
||||
@ -1620,8 +1618,7 @@ class SumkDFT(object):
|
||||
for ik in range(self.n_k):
|
||||
n_orb = self.n_orbitals[ik, ind]
|
||||
MMat = np.identity(n_orb, complex)
|
||||
MMat = self.hopping[
|
||||
ik, ind, 0:n_orb, 0:n_orb] - (1 - 2 * isp) * self.h_field * MMat
|
||||
MMat = self.hopping[ik, ind, 0:n_orb, 0:n_orb] - (1 - 2 * isp) * self.h_field * MMat
|
||||
projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb]
|
||||
self.Hsumk[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat),
|
||||
projmat.conjugate().transpose())
|
||||
@ -2162,52 +2159,56 @@ class SumkDFT(object):
|
||||
if dm_type is None:
|
||||
dm_type = self.dft_code
|
||||
|
||||
assert dm_type in ('vasp', 'wien2k','elk', 'qe'), "'dm_type' must be either 'vasp', 'wienk', 'elk' or 'qe'"
|
||||
#default file names
|
||||
assert dm_type in ('vasp', 'wien2k', 'elk', 'qe'), "'dm_type' must be either 'vasp', 'wienk', 'elk' or 'qe'"
|
||||
# default file names
|
||||
if filename is None:
|
||||
if dm_type == 'wien2k':
|
||||
filename = 'dens_mat.dat'
|
||||
elif dm_type == 'vasp':
|
||||
filename = 'GAMMA'
|
||||
# use new h5 interface to vasp by default, if not wanted specify dm_type='vasp' + filename='GAMMA'
|
||||
filename = 'vaspgamma.h5'
|
||||
elif dm_type == 'elk':
|
||||
filename = 'DMATDMFT.OUT'
|
||||
elif dm_type == 'qe':
|
||||
filename = self.hdf_file
|
||||
|
||||
|
||||
assert isinstance(filename, str), ("calc_density_correction: "
|
||||
"filename has to be a string!")
|
||||
"filename has to be a string!")
|
||||
|
||||
assert kpts_to_write is None or dm_type == 'vasp', ('Selecting k-points only'
|
||||
+'implemented for vasp')
|
||||
+ 'implemented for vasp')
|
||||
|
||||
ntoi = self.spin_names_to_ind[self.SO]
|
||||
spn = self.spin_block_names[self.SO]
|
||||
dens = {sp: 0.0 for sp in spn}
|
||||
band_en_correction = 0.0
|
||||
|
||||
# Fetch Fermi weights and energy window band indices
|
||||
if dm_type in ['vasp','qe']:
|
||||
# Fetch Fermi weights and energy window band indices
|
||||
if dm_type in ['vasp', 'qe']:
|
||||
fermi_weights = 0
|
||||
band_window = 0
|
||||
n_k_ibz = self.n_k
|
||||
if mpi.is_master_node():
|
||||
with HDFArchive(self.hdf_file,'r') as ar:
|
||||
with HDFArchive(self.hdf_file, 'r') as ar:
|
||||
fermi_weights = ar['dft_misc_input']['dft_fermi_weights']
|
||||
band_window = ar['dft_misc_input']['band_window']
|
||||
|
||||
if 'n_k_ibz' in ar['dft_misc_input']:
|
||||
n_k_ibz = ar['dft_misc_input']['n_k_ibz']
|
||||
fermi_weights = mpi.bcast(fermi_weights)
|
||||
band_window = mpi.bcast(band_window)
|
||||
n_k_ibz = mpi.bcast(n_k_ibz)
|
||||
|
||||
# Convert Fermi weights to a density matrix
|
||||
# Convert Fermi weights to a density matrix
|
||||
dens_mat_dft = {}
|
||||
for sp in spn:
|
||||
dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(complex) for ik in range(self.n_k)]
|
||||
|
||||
|
||||
# Set up deltaN:
|
||||
deltaN = {}
|
||||
for sp in spn:
|
||||
deltaN[sp] = [np.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[
|
||||
ik, ntoi[sp]]], complex) for ik in range(self.n_k)]
|
||||
ik, ntoi[sp]]], complex) for ik in range(self.n_k)]
|
||||
|
||||
ikarray = np.arange(self.n_k)
|
||||
for ik in mpi.slice_array(ikarray):
|
||||
@ -2218,7 +2219,7 @@ class SumkDFT(object):
|
||||
for bname, gf in G_latt:
|
||||
G_latt_rot = gf.copy()
|
||||
G_latt_rot << self.upfold(
|
||||
ik, 0, bname, G_latt[bname], gf,shells='csc')
|
||||
ik, 0, bname, G_latt[bname], gf, shells='csc')
|
||||
|
||||
G_latt[bname] = G_latt_rot.copy()
|
||||
|
||||
@ -2229,16 +2230,16 @@ class SumkDFT(object):
|
||||
dens[bname] += self.bz_weights[ik] * G_latt[bname].total_density()
|
||||
else:
|
||||
dens[bname] += self.bz_weights[ik] * G_latt[bname].total_density(beta)
|
||||
if dm_type in ['vasp','qe']:
|
||||
# In 'vasp'-mode subtract the DFT density matrix
|
||||
if dm_type in ['vasp', 'qe']:
|
||||
# In 'vasp'-mode subtract the DFT density matrix
|
||||
nb = self.n_orbitals[ik, ntoi[bname]]
|
||||
diag_inds = np.diag_indices(nb)
|
||||
deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik][:nb]
|
||||
|
||||
if self.charge_mixing and self.deltaNOld is not None:
|
||||
G2 = np.sum(self.kpts_cart[ik,:]**2)
|
||||
G2 = np.sum(self.kpts_cart[ik, :] ** 2)
|
||||
# Kerker mixing
|
||||
mix_fac = self.charge_mixing_alpha * G2 / (G2 + self.charge_mixing_gamma**2)
|
||||
mix_fac = self.charge_mixing_alpha * G2 / (G2 + self.charge_mixing_gamma ** 2)
|
||||
deltaN[bname][ik][diag_inds] = (1.0 - mix_fac) * self.deltaNOld[bname][ik][diag_inds] + mix_fac * deltaN[bname][ik][diag_inds]
|
||||
dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real
|
||||
isp = ntoi[bname]
|
||||
@ -2283,10 +2284,8 @@ class SumkDFT(object):
|
||||
f.write("%s\n" % self.n_orbitals[ik, 0])
|
||||
for inu in range(self.n_orbitals[ik, 0]):
|
||||
for imu in range(self.n_orbitals[ik, 0]):
|
||||
valre = (deltaN['up'][ik][
|
||||
inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
|
||||
valim = (deltaN['up'][ik][
|
||||
inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
|
||||
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
|
||||
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
|
||||
f.write("%.14f %.14f " % (valre, valim))
|
||||
f.write("\n")
|
||||
f.write("\n")
|
||||
@ -2305,32 +2304,53 @@ class SumkDFT(object):
|
||||
fout.write("%s\n" % self.n_orbitals[ik, isp])
|
||||
for inu in range(self.n_orbitals[ik, isp]):
|
||||
for imu in range(self.n_orbitals[ik, isp]):
|
||||
fout.write("%.14f %.14f " % (deltaN[sp][ik][
|
||||
inu, imu].real, deltaN[sp][ik][inu, imu].imag))
|
||||
fout.write("%.14f %.14f " % (deltaN[sp][ik][inu, imu].real, deltaN[sp][ik][inu, imu].imag))
|
||||
fout.write("\n")
|
||||
fout.write("\n")
|
||||
fout.close()
|
||||
elif dm_type == 'vasp':
|
||||
if kpts_to_write is None:
|
||||
if kpts_to_write is None and self.n_k == n_k_ibz:
|
||||
kpts_to_write = np.arange(self.n_k)
|
||||
elif kpts_to_write is None and n_k_ibz < self.n_k:
|
||||
# If the number of IBZ k-points is less than the total number of k-points,
|
||||
# then we only write the IBZ k-points, which are the first n_k_ibz k-points
|
||||
# in VASP
|
||||
kpts_to_write = np.arange(n_k_ibz)
|
||||
else:
|
||||
assert np.min(kpts_to_write) >= 0 and np.max(kpts_to_write) < self.n_k
|
||||
|
||||
assert self.SP == 0, "Spin-polarized density matrix is not implemented"
|
||||
|
||||
if mpi.is_master_node():
|
||||
with open(filename, 'w') as f:
|
||||
f.write(" %i -1 ! Number of k-points, default number of bands\n"%len(kpts_to_write))
|
||||
for index, ik in enumerate(kpts_to_write):
|
||||
ib1 = band_window[0][ik, 0]
|
||||
ib2 = band_window[0][ik, 1]
|
||||
f.write(" %i %i %i\n"%(index + 1, ib1, ib2))
|
||||
for inu in range(self.n_orbitals[ik, 0]):
|
||||
for imu in range(self.n_orbitals[ik, 0]):
|
||||
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
|
||||
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
|
||||
f.write(" %.14f %.14f"%(valre, valim))
|
||||
f.write("\n")
|
||||
if filename == 'vaspgamma.h5':
|
||||
with HDFArchive('vaspgamma.h5', 'w') as vasp_h5:
|
||||
# only store the ibz kpoints in the h5
|
||||
bnd_win_towrite = [band_window[0][:n_k_ibz,:]]
|
||||
vasp_h5['band_window'] = bnd_win_towrite
|
||||
vasp_h5.create_group('deltaN')
|
||||
vasp_h5['deltaN']['up'] = deltaN['up'][:n_k_ibz]
|
||||
vasp_h5['deltaN']['down'] = deltaN['down'][:n_k_ibz]
|
||||
else:
|
||||
with open(filename, 'w') as f:
|
||||
f.write(" -1 -1 ! Number of k-points, default number of bands\n") # % len(kpts_to_write))
|
||||
for index, ik in enumerate(kpts_to_write):
|
||||
ib1 = band_window[0][ik, 0]
|
||||
ib2 = band_window[0][ik, 1]
|
||||
f.write(" %i %i %i\n" % (index + 1, ib1, ib2))
|
||||
for inu in range(self.n_orbitals[ik, 0]):
|
||||
for imu in range(self.n_orbitals[ik, 0]):
|
||||
if (self.SO == 1):
|
||||
valre = (deltaN['ud'][ik][inu, imu].real) / 1.0
|
||||
valim = (deltaN['ud'][ik][inu, imu].imag) / 1.0
|
||||
f.write(" %.14f %.14f" % (valre, valim))
|
||||
else:
|
||||
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][
|
||||
inu, imu].real) / 2.0
|
||||
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][
|
||||
inu, imu].imag) / 2.0
|
||||
f.write(" %.14f %.14f" % (valre, valim))
|
||||
f.write("\n")
|
||||
|
||||
|
||||
elif dm_type == 'elk':
|
||||
# output each k-point density matrix for Elk
|
||||
@ -2437,8 +2457,7 @@ class SumkDFT(object):
|
||||
dim = self.corr_shells[icrsh]['dim']
|
||||
n_orb = self.n_orbitals[ik, 0]
|
||||
projmat = self.proj_mat[ik, 0, icrsh, 0:dim, 0:n_orb]
|
||||
dens_mat[icrsh][
|
||||
:, :] += np.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
|
||||
dens_mat[icrsh][:, :] += np.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
|
||||
|
||||
if self.symm_op != 0:
|
||||
dens_mat = self.symmcorr.symmetrize(dens_mat)
|
||||
@ -2448,8 +2467,7 @@ class SumkDFT(object):
|
||||
for icrsh in range(self.n_corr_shells):
|
||||
if self.rot_mat_time_inv[icrsh] == 1:
|
||||
dens_mat[icrsh] = dens_mat[icrsh].conjugate()
|
||||
dens_mat[icrsh] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]),
|
||||
self.rot_mat[icrsh])
|
||||
dens_mat[icrsh] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]),self.rot_mat[icrsh])
|
||||
|
||||
return dens_mat
|
||||
|
||||
|
Binary file not shown.
Binary file not shown.
Binary file not shown.
16
test/python/plovasp/converter/svo.cfg
Normal file
16
test/python/plovasp/converter/svo.cfg
Normal file
@ -0,0 +1,16 @@
|
||||
[General]
|
||||
BASENAME = converter/svo
|
||||
|
||||
[Group 1]
|
||||
SHELLS = 1
|
||||
NORMALIZE = True
|
||||
EWINDOW = -1.4 2.0
|
||||
|
||||
[Shell 1]
|
||||
LSHELL = 2
|
||||
IONS = 2
|
||||
|
||||
TRANSFORM = 1.0 0.0 0.0 0.0 0.0
|
||||
0.0 1.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 1.0 0.0
|
||||
|
BIN
test/python/plovasp/converter/svo.ref.h5
Normal file
BIN
test/python/plovasp/converter/svo.ref.h5
Normal file
Binary file not shown.
BIN
test/python/plovasp/converter/svo/vaspout.h5
Normal file
BIN
test/python/plovasp/converter/svo/vaspout.h5
Normal file
Binary file not shown.
42
test/python/plovasp/converter/test_converter_svo.py
Normal file
42
test/python/plovasp/converter/test_converter_svo.py
Normal file
@ -0,0 +1,42 @@
|
||||
|
||||
import os
|
||||
import rpath
|
||||
_rpath = os.path.dirname(rpath.__file__) + '/'
|
||||
|
||||
from triqs_dft_tools.converters.plovasp.converter import generate_and_output_as_text
|
||||
from triqs_dft_tools.converters import VaspConverter
|
||||
import mytest
|
||||
|
||||
################################################################################
|
||||
#
|
||||
# TestConverterOneSite
|
||||
#
|
||||
################################################################################
|
||||
class TestConverterSVO(mytest.MyTestCase):
|
||||
"""
|
||||
Function:
|
||||
|
||||
def generate_and_output_as_text(pars, el_struct)
|
||||
and
|
||||
VaspConverter
|
||||
|
||||
Scenarios:
|
||||
|
||||
- Parse config file and produce a correct converted h5-file
|
||||
"""
|
||||
# Scenario 1
|
||||
def test_convert_svo(self):
|
||||
generate_and_output_as_text(_rpath + 'svo.cfg', _rpath + 'svo/')
|
||||
|
||||
test_file = _rpath + 'svo.test.h5'
|
||||
converter = VaspConverter(filename=_rpath + 'svo',
|
||||
hdf_filename=test_file)
|
||||
|
||||
converter.convert_dft_input()
|
||||
|
||||
expected_file = _rpath + 'svo.ref.h5'
|
||||
self.assertH5FileEqual(test_file, expected_file)
|
||||
|
||||
if __name__ == '__main__':
|
||||
import unittest
|
||||
unittest.main(verbosity=2, buffer=False)
|
Loading…
x
Reference in New Issue
Block a user