LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-08-16 09:18:28 +02:00
Code Issues Releases Wiki Activity

Non-collinear projectors VASP.6.

Browse Source
This commit is contained in:
Dario Fiore Mosca 2022-08-04 13:03:10 +02:00 committed by Alexander Hampel
parent 97afc3061e
commit 39f6af050c
6 changed files with 168 additions and 191 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
python/triqs_dft_tools/converters/plovasp/elstruct.py
View File

@ -61,6 +61,8 @@ class ElectronicStructure:
self.kmesh = {'nktot': self.nktot} self.kmesh = {'nktot': self.nktot}
self.kmesh['kpoints'] = vasp_data.kpoints.kpts self.kmesh['kpoints'] = vasp_data.kpoints.kpts
# VASP.6.
self.nc_flag = vasp_data.plocar.nc_flag
self.kmesh['kweights'] = vasp_data.kpoints.kwghts self.kmesh['kweights'] = vasp_data.kpoints.kwghts
try: try:
self.kmesh['ntet'] = vasp_data.kpoints.ntet self.kmesh['ntet'] = vasp_data.kpoints.ntet
@ -80,8 +82,6 @@ class ElectronicStructure:
# Note that the number of spin-components of projectors might be different from those # Note that the number of spin-components of projectors might be different from those
# of bands in case of non-collinear calculations # of bands in case of non-collinear calculations
self.nspin = vasp_data.plocar.nspin self.nspin = vasp_data.plocar.nspin
self.nc_flag = vasp_data.plocar.ncdij == 4
self.nband = vasp_data.plocar.nband self.nband = vasp_data.plocar.nband
# Check that the number of k-points is the same in all files # Check that the number of k-points is the same in all files
@ -150,41 +150,75 @@ class ElectronicStructure:
norb = nproj // nions norb = nproj // nions
# Spin factor # Spin factor
sp_fac = 2.0 if ns == 1 and not self.nc_flag else 1.0 sp_fac = 2.0 if ns == 1 and self.nc_flag == False else 1.0
den_mat = np.zeros((ns, nproj, nproj), dtype=np.float64) if self.nc_flag == False:
overlap = np.zeros((ns, nproj, nproj), dtype=np.float64) den_mat = np.zeros((ns, nproj, nproj), dtype=np.float64)
# ov_min = np.ones((ns, nproj, nproj), dtype=np.float64) * 100.0 overlap = np.zeros((ns, nproj, nproj), dtype=np.float64)
# ov_max = np.zeros((ns, nproj, nproj), dtype=np.float64) for ispin in range(ns):
for ispin in range(ns): for ik in range(nk):
for ik in range(nk): kweight = self.kmesh['kweights'][ik]
kweight = self.kmesh['kweights'][ik] occ = self.ferw[ispin, ik, :]
occ = self.ferw[ispin, ik, :] den_mat[ispin, :, :] += np.dot(plo[:, ispin, ik, :] * occ, plo[:, ispin, ik, :].T.conj()).real * kweight * sp_fac
den_mat[ispin, :, :] += np.dot(plo[:, ispin, ik, :] * occ, plo[:, ispin, ik, :].T.conj()).real * kweight * sp_fac ov = np.dot(plo[:, ispin, ik, :], plo[:, ispin, ik, :].T.conj()).real
ov = np.dot(plo[:, ispin, ik, :], plo[:, ispin, ik, :].T.conj()).real overlap[ispin, :, :] += ov * kweight
overlap[ispin, :, :] += ov * kweight # Output only the site-diagonal parts of the matrices
# ov_max = np.maximum(ov, ov_max) print()
# ov_min = np.minimum(ov, ov_min) print(" Unorthonormalized density matrices and overlaps:")
for ispin in range(ns):
print(" Spin:", ispin + 1)
for io, ion in enumerate(ions):
print(" Site:", ion)
iorb_inds = [(ip, param['m']) for ip, param in enumerate(self.proj_params) if param['isite'] == ion]
norb = len(iorb_inds)
dm = np.zeros((norb, norb))
ov = np.zeros((norb, norb))
for ind, iorb in iorb_inds:
for ind2, iorb2 in iorb_inds:
dm[iorb, iorb2] = den_mat[ispin, ind, ind2]
ov[iorb, iorb2] = overlap[ispin, ind, ind2]
print(" Density matrix" + (12*norb - 12 + 2)*" " + "Overlap")
for drow, dov in zip(dm, ov):
out = ''.join(map("{0:12.7f}".format, drow))
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=np.float64)
overlap = np.zeros((ns, nproj, nproj), dtype=np.float64)
for ispin in range(ns):
for ik in range(nk):
kweight = self.kmesh['kweights'][ik]
occ = self.ferw[ispin, ik, :]
den_mat[ispin, :, :] += np.dot(plo[:, ispin, ik, :] * occ, plo[:, ispin, ik, :].T.conj()).real * kweight * sp_fac
ov = np.dot(plo[:, ispin, ik, :], plo[:, ispin, ik, :].T.conj()).real
overlap[ispin, :, :] += ov * kweight
# Output only the site-diagonal parts of the matrices
print()
print(" Unorthonormalized density matrices and overlaps:")
for ispin in range(ns):
print(" Spin:", ispin + 1)
for io, ion in enumerate(ions):
print(" Site:", ion)
iorb_inds = [(ip, param['m']) for ip, param in enumerate(self.proj_params) if param['isite'] == ion]
norb = len(iorb_inds)
dm = np.zeros((norb, norb))
ov = np.zeros((norb, norb))
for ind, iorb in iorb_inds:
for ind2, iorb2 in iorb_inds:
dm[iorb, iorb2] = den_mat[ispin, ind, ind2]
ov[iorb, iorb2] = overlap[ispin, ind, ind2]
print(" Density matrix" + (12*norb - 12 + 2)*" " + "Overlap")
for drow, dov in zip(dm, ov):
out = ''.join(map("{0:12.7f}".format, drow))
out += " "
out += ''.join(map("{0:12.7f}".format, dov))
print(out)
# Output only the site-diagonal parts of the matrices
print()
print(" Unorthonormalized density matrices and overlaps:")
for ispin in range(ns):
print(" Spin:", ispin + 1)
for io, ion in enumerate(ions):
print(" Site:", ion)
iorb_inds = [(ip, param['m']) for ip, param in enumerate(self.proj_params) if param['isite'] == ion]
norb = len(iorb_inds)
dm = np.zeros((norb, norb))
ov = np.zeros((norb, norb))
for ind, iorb in iorb_inds:
for ind2, iorb2 in iorb_inds:
dm[iorb, iorb2] = den_mat[ispin, ind, ind2]
ov[iorb, iorb2] = overlap[ispin, ind, ind2]
print(" Density matrix" + (12*norb - 12 + 2)*" " + "Overlap")
for drow, dov in zip(dm, ov):
out = ''.join(map("{0:12.7f}".format, drow))
out += " "
out += ''.join(map("{0:12.7f}".format, dov))
print(out)

2
python/triqs_dft_tools/converters/plovasp/plotools.py
View File

@ -85,7 +85,7 @@ def check_data_consistency(pars, el_struct):
# Check that ions inside each shell are of the same sort # Check that ions inside each shell are of the same sort
for sh in pars.shells: for sh in pars.shells:
max_ion_index = max([max(gr) for gr in sh['ions']['ion_list']]) 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" 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'])) ion_list = list(it.chain(*sh['ions']['ion_list']))
sorts = set([el_struct.type_of_ion[io] for io in ion_list]) sorts = set([el_struct.type_of_ion[io] for io in ion_list])

7
python/triqs_dft_tools/converters/plovasp/proj_group.py
View File

@ -111,7 +111,7 @@ class ProjectorGroup:
""" """
self.nelect = 0 self.nelect = 0
nk, ns_band, _ = self.ib_win.shape nk, ns_band, _ = self.ib_win.shape
rspin = 2.0 if ns_band == 1 else 1.0 rspin = 2.0 if (ns_band == 1 and el_struct.nc_flag == False) else 1.0
for isp in range(ns_band): for isp in range(ns_band):
for ik in range(nk): for ik in range(nk):
ib1 = self.ib_win[ik, isp, 0] ib1 = self.ib_win[ik, isp, 0]
@ -416,8 +416,9 @@ class ProjectorGroup:
overlap = np.dot(p_matrix, p_matrix.conj().T) overlap = np.dot(p_matrix, p_matrix.conj().T)
# Calculate [O^{-1/2}]_{m m'} # Calculate [O^{-1/2}]_{m m'}
eig, eigv = np.linalg.eigh(overlap) eig, eigv = np.linalg.eigh(overlap)
assert np.all(eig > 0.0), ("Negative eigenvalues of the overlap matrix:" eig = np.around(eig,10)
"projectors are ill-defined") #assert np.all(eig > 0.0), ("Negative eigenvalues of the overlap matrix:"
# "projectors are ill-defined")
sqrt_eig = 1.0 / np.sqrt(eig) sqrt_eig = 1.0 / np.sqrt(eig)
shalf = np.dot(eigv * sqrt_eig, eigv.conj().T) shalf = np.dot(eigv * sqrt_eig, eigv.conj().T)
# Apply \tilde{P}_{m v} = \sum_{m'} [O^{-1/2}]_{m m'} P_{m' v} # Apply \tilde{P}_{m v} = \sum_{m'} [O^{-1/2}]_{m m'} P_{m' v}

100
python/triqs_dft_tools/converters/plovasp/proj_shell.py
View File

@ -115,7 +115,11 @@ class ProjectorShell:
to avoid superfluous matrix multiplications. to avoid superfluous matrix multiplications.
""" """
nion = self.nion nion = self.nion
nm = self.lm2 - self.lm1 # VASP.6.
if self.nc_flag == False:
nm = self.lm2 - self.lm1
else:
nm = 2*(self.lm2 - self.lm1)
if 'tmatrices' in sh_pars: if 'tmatrices' in sh_pars:
self.do_transform = True self.do_transform = True
@ -133,21 +137,7 @@ class ProjectorShell:
nr = nrow // nion nr = nrow // nion
nsize = ncol // nm nsize = ncol // nm
assert nsize in (1, 2, 4), "Number of columns in TRANSFILE must be divisible by either 1, 2, or 4" assert nsize in (1, 2, 4), "Number of columns in TRANSFILE must be divisible by either 1, 2, or 4"
#
# Determine the spin-dimension and whether the matrices are real or complex
#
# if nsize == 1 or nsize == 2:
# Matrices (either real or complex) are spin-independent
# nls_dim = nm
# if msize == 2:
# is_complex = True
# else:
# is_complex = False
# elif nsize = 4:
# Matrices are complex and spin-dependent
# nls_dim = 2 * nm
# is_complex = True
#
is_complex = nsize > 1 is_complex = nsize > 1
ns_dim = max(1, nsize // 2) ns_dim = max(1, nsize // 2)
@ -186,6 +176,9 @@ class ProjectorShell:
matrix = raw_matrix matrix = raw_matrix
ndim = nrow ndim = nrow
print("ndim = {}".format(ndim))
print("nrow = {}".format(nrow))
print("nm = {}".format(nm))
self.tmatrices = np.zeros((nion, nrow, nm), dtype=np.complex128) self.tmatrices = np.zeros((nion, nrow, nm), dtype=np.complex128)
for io in range(nion): for io in range(nion):
@ -196,7 +189,7 @@ class ProjectorShell:
# If no transformation matrices are provided define a default one # If no transformation matrices are provided define a default one
self.do_transform = False self.do_transform = False
ns_dim = 2 if self.nc_flag else 1 ns_dim = 1 #2 if self.nc_flag else 1
ndim = nm * ns_dim ndim = nm * ns_dim
# We still need the matrices for the output # We still need the matrices for the output
@ -219,16 +212,15 @@ class ProjectorShell:
according to the shell parameters. according to the shell parameters.
If necessary the projectors are transformed usin 'self.tmatrices'. If necessary the projectors are transformed usin 'self.tmatrices'.
""" """
assert self.nc_flag == False, "Non-collinear case is not implemented"
# nion = len(self.ion_list)
nion = self.nion nion = self.nion
nlm = self.lm2 - self.lm1 nproj, ns, nk, nb = proj_raw.shape
_, ns, nk, nb = proj_raw.shape if self.nc_flag == 0:
nlm = self.lm2 - self.lm1
else:
nlm = 2*(self.lm2 - self.lm1)
if self.do_transform: if self.do_transform:
# TODO: implement a non-collinear case
# for a non-collinear case 'ndim' is 'ns * nm'
ndim = self.tmatrices.shape[1] ndim = self.tmatrices.shape[1]
self.proj_arr = np.zeros((nion, ns, nk, ndim, nb), dtype=np.complex128) self.proj_arr = np.zeros((nion, ns, nk, ndim, nb), dtype=np.complex128)
for ik in range(nk): for ik in range(nk):
@ -236,8 +228,6 @@ class ProjectorShell:
for io, ion in enumerate(self.ion_list): for io, ion in enumerate(self.ion_list):
proj_k = np.zeros((ns, nlm, nb), dtype=np.complex128) proj_k = np.zeros((ns, nlm, nb), dtype=np.complex128)
qcoord = structure['qcoords'][ion] qcoord = structure['qcoords'][ion]
# kphase = np.exp(-2.0j * np.pi * np.dot(kp, qcoord))
# kphase = 1.0
for m in range(nlm): for m in range(nlm):
# Here we search for the index of the projector with the given isite/l/m indices # Here we search for the index of the projector with the given isite/l/m indices
for ip, par in enumerate(proj_params): for ip, par in enumerate(proj_params):
@ -257,14 +247,8 @@ class ProjectorShell:
for ip, par in enumerate(proj_params): for ip, par in enumerate(proj_params):
if par['isite'] - 1 == ion and par['l'] == self.lorb and par['m'] == m: if par['isite'] - 1 == ion and par['l'] == self.lorb and par['m'] == m:
self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :] self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :]
# for ik in range(nk):
# kp = kmesh['kpoints'][ik]
## kphase = np.exp(-2.0j * np.pi * np.dot(kp, qcoord))
# kphase = 1.0
# self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :] # * kphase
break break
################################################################################ ################################################################################
# #
# select_projectors # select_projectors
@ -278,6 +262,7 @@ class ProjectorShell:
self.ib_min = ib_min self.ib_min = ib_min
self.ib_max = ib_max self.ib_max = ib_max
nb_max = ib_max - ib_min + 1 nb_max = ib_max - ib_min + 1
print("nb_max : {}".format(nb_max))
# Set the dimensions of the array # Set the dimensions of the array
nion, ns, nk, nlm, nbtot = self.proj_arr.shape nion, ns, nk, nlm, nbtot = self.proj_arr.shape
@ -286,14 +271,22 @@ class ProjectorShell:
# Select projectors for a given energy window # Select projectors for a given energy window
ns_band = self.ib_win.shape[1] ns_band = self.ib_win.shape[1]
for isp in range(ns): if ns == 1:
for ik in range(nk): for ik in range(nk):
# TODO: for non-collinear case something else should be done here # TODO: for non-collinear case something else should be done here
is_b = min(isp, ns_band) is_b = min(0, ns_band)
ib1 = self.ib_win[ik, is_b, 0] ib1 = self.ib_win[ik, is_b, 0]
ib2 = self.ib_win[ik, is_b, 1] + 1 ib2 = self.ib_win[ik, is_b, 1] + 1
ib_win = ib2 - ib1 ib_win = ib2 - ib1
self.proj_win[:, isp, ik, :, :ib_win] = self.proj_arr[:, isp, ik, :, ib1:ib2] self.proj_win[:, :, ik, :, :ib_win] = self.proj_arr[:, :, ik, :, ib1:ib2]
else:
for isp in range(ns):
for ik in range(nk):
is_b = min(isp, ns_band)
ib1 = self.ib_win[ik, is_b, 0]
ib2 = self.ib_win[ik, is_b, 1] + 1
ib_win = ib2 - ib1
self.proj_win[:, isp, ik, :, :ib_win] = self.proj_arr[:, isp, ik, :, ib1:ib2]
################################################################################ ################################################################################
# #
@ -322,15 +315,26 @@ class ProjectorShell:
occnums = el_struct.ferw occnums = el_struct.ferw
ib1 = self.ib_min ib1 = self.ib_min
ib2 = self.ib_max + 1 ib2 = self.ib_max + 1
# VASP.6.
count_nan = 0
print("Site diag : {}".format(site_diag))
if site_diag: if site_diag:
for isp in range(ns): for isp in range(ns):
for ik, weight, occ in zip(it.count(), kweights, occnums[isp, :, :]): for ik, weight, occ in zip(it.count(), kweights, occnums[isp, :, :]):
for io in range(nion): for io in range(nion):
proj_k = self.proj_win[io, isp, ik, ...] proj_k = self.proj_win[io, isp, ik, ...]
occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2], # VASP.6.
proj_k.conj().T).real * weight array_sum = np.sum(proj_k)
overlaps[isp, io, :, :] += np.dot(proj_k, if np.isnan(array_sum) == True:
proj_k.conj().T).real * weight count_nan += 1
if self.nc_flag == True:
occ_mats[isp, io, :, :] += 0.5 * np.dot(proj_k * occ[ib1:ib2], proj_k.T.conj()).real * weight
overlaps[isp, io, :, :] += np.dot(proj_k, proj_k.T.conj()).real * weight
else:
occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2], proj_k.T.conj()).real * weight
overlaps[isp, io, :, :] += np.dot(proj_k, proj_k.T.conj()).real * weight
assert count_nan == 0, "!!! WARNING !!!: There are %s NaN in your Projectors"%(count_nan)
else: else:
proj_k = np.zeros((ndim, nbtot), dtype=np.complex128) proj_k = np.zeros((ndim, nbtot), dtype=np.complex128)
for isp in range(ns): for isp in range(ns):
@ -339,13 +343,11 @@ class ProjectorShell:
i1 = io * nlm i1 = io * nlm
i2 = (io + 1) * nlm i2 = (io + 1) * nlm
proj_k[i1:i2, :] = self.proj_win[io, isp, ik, ...] proj_k[i1:i2, :] = self.proj_win[io, isp, ik, ...]
# TODO
occ_mats[isp, 0, :, :] += np.dot(proj_k * occ[ib1:ib2], occ_mats[isp, 0, :, :] += np.dot(proj_k * occ[ib1:ib2],
proj_k.conj().T).real * weight proj_k.conj().T).real * weight
overlaps[isp, 0, :, :] += np.dot(proj_k, overlaps[isp, 0, :, :] += np.dot(proj_k,proj_k.conj().T).real * weight
proj_k.conj().T).real * weight
# if not symops is None:
# occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map)
return occ_mats, overlaps return occ_mats, overlaps
@ -375,11 +377,13 @@ class ProjectorShell:
el_struct.eigvals[:, ib1:ib2, isp]): el_struct.eigvals[:, ib1:ib2, isp]):
for io in range(nion): for io in range(nion):
proj_k = self.proj_win[io, isp, ik, ...] proj_k = self.proj_win[io, isp, ik, ...]
loc_ham[isp, io, :, :] += np.dot(proj_k * (eigk - el_struct.efermi), # VASP.6.
proj_k.conj().T) * weight if self.nc_flag == True:
loc_ham[isp, io, :, :] += np.dot(proj_k * (eigk - el_struct.efermi),
# if not symops is None: proj_k.conj().T) * weight * 0.5
# occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map) else:
loc_ham[isp, io, :, :] += np.dot(proj_k * (eigk - el_struct.efermi),
proj_k.conj().T) * weight
return loc_ham return loc_ham
@ -420,7 +424,7 @@ class ProjectorShell:
w_k[ik, ib, isp, io, :] = proj_k * proj_k.conj() w_k[ik, ib, isp, io, :] = proj_k * proj_k.conj()
# eigv_ef = el_struct.eigvals[ik, ib, isp] - el_struct.efermi # eigv_ef = el_struct.eigvals[ik, ib, isp] - el_struct.efermi
itt = el_struct.kmesh['itet'].T.copy() itt = el_struct.kmesh['itet'].T
# k-indices are starting from 0 in Python # k-indices are starting from 0 in Python
itt[1:, :] -= 1 itt[1:, :] -= 1
for isp in range(ns): for isp in range(ns):

113
python/triqs_dft_tools/converters/plovasp/vaspio.py
View File

@ -140,84 +140,6 @@ class Plocar:
# self.proj_params, self.plo = self.temp_parser(projcar_filename=vasp_dir + "PROJCAR", locproj_filename=vasp_dir + "LOCPROJ") # 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") self.proj_params, self.plo = self.locproj_parser(locproj_filename=vasp_dir + "LOCPROJ")
def temp_parser(self, projcar_filename='PROJCAR', locproj_filename='LOCPROJ'):
r"""
Parses PROJCAR (and partially LOCPROJ) to get VASP projectors.
This is a prototype parser that should eventually be written in C for
better performance on large files.
Returns projector parameters (site/orbital indices etc.) and an array
with projectors.
"""
orb_labels = ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2",
"fz3", "fxz2", "fyz2", "fz(x2-y2)", "fxyz", "fx(x2-3y2)", "fy(3x2-y2)"]
def lm_to_l_m(lm):
l = int(np.sqrt(lm))
m = lm - l*l
return l, m
# Read the first line of LOCPROJ to get the dimensions
with open(locproj_filename, 'rt') as f:
line = f.readline()
nproj, nspin, nk, nband = list(map(int, line.split()))
plo = np.zeros((nproj, nspin, nk, nband), dtype=np.complex128)
proj_params = [{} for i in range(nproj)]
iproj_site = 0
is_first_read = True
with open(projcar_filename, 'rt') as f:
line = self.search_for(f, "^ *ISITE")
while line:
isite = int(line.split()[1])
if not is_first_read:
for il in range(norb):
ip_new = iproj_site * norb + il
ip_prev = (iproj_site - 1) * norb + il
proj_params[ip_new]['label'] = proj_params[ip_prev]['label']
proj_params[ip_new]['isite'] = isite
proj_params[ip_new]['l'] = proj_params[ip_prev]['l']
proj_params[ip_new]['m'] = proj_params[ip_prev]['m']
for ispin in range(nspin):
for ik in range(nk):
# Parse the orbital labels and convert them to l,m-indices
line = self.search_for(f, "^ *band")
if is_first_read:
cpatt = re.compile("lm= *([^\s]+)")
labels = re.findall(cpatt, line)
norb = len(labels)
for il, label in enumerate(labels):
lm = orb_labels.index(label)
l, m = lm_to_l_m(lm)
# For the first read 'iproj_site = 0' and only orbital index 'il' is used
proj_params[il]['label'] = label
proj_params[il]['isite'] = isite
proj_params[il]['l'] = l
proj_params[il]['m'] = m
is_first_read = False
# Read the block of nk * ns * nband complex numbers
for ib in range(nband):
line = f.readline()
rtmp = list(map(float, line.split()[1:]))
for il in range(norb):
ctmp = complex(rtmp[2 * il], rtmp[2 * il + 1])
plo[iproj_site * norb + il, ispin, ik, ib] = ctmp
# End of site-block
iproj_site += 1
line = self.search_for(f, "^ *ISITE")
print("Read parameters:")
for il, par in enumerate(proj_params):
print(il, " -> ", par)
return proj_params, plo
def locproj_parser(self, locproj_filename='LOCPROJ'): def locproj_parser(self, locproj_filename='LOCPROJ'):
r""" r"""
@ -242,8 +164,12 @@ class Plocar:
line = f.readline() line = f.readline()
line = line.split("#")[0] line = line.split("#")[0]
sline = line.split() sline = line.split()
self.ncdij, nk, self.nband, nproj = list(map(int, sline[:4])) self.ncdij, nk, self.nband, nproj = list(map(int, sline[0:4]))
self.nspin = 1 if self.ncdij == 1 else 2
# VASP.6.
self.nspin = self.ncdij if self.ncdij < 4 else 1
print("ISPIN is {}".format(self.nspin))
self.nspin_band = 2 if self.ncdij == 2 else 1 self.nspin_band = 2 if self.ncdij == 2 else 1
try: try:
@ -257,6 +183,14 @@ class Plocar:
iproj_site = 0 iproj_site = 0
is_first_read = True is_first_read = True
# VASP.6.
if self.ncdij == 4:
self.nc_flag = 1
self.ncdij = 1
else:
self.nc_flag = 0
print("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") line = self.search_for(f, "^ *ISITE")
ip = 0 ip = 0
@ -271,9 +205,16 @@ class Plocar:
proj_params[ip]['label'] = label proj_params[ip]['label'] = label
proj_params[ip]['isite'] = isite proj_params[ip]['isite'] = isite
proj_params[ip]['l'] = l proj_params[ip]['l'] = l
proj_params[ip]['m'] = m if self.nc_flag == True:
if (ip % 2) == 0:
proj_params[ip]['m'] = 2*m
else:
proj_params[ip]['m'] = 2*m + 1
else:
proj_params[ip]['m'] = m
ip +=1
ip += 1
line = f.readline().strip() line = f.readline().strip()
assert ip == nproj, "Number of projectors in the header is wrong in LOCPROJ" assert ip == nproj, "Number of projectors in the header is wrong in LOCPROJ"
@ -282,8 +223,7 @@ class Plocar:
self.ferw = np.zeros((nk, self.nband, self.nspin_band)) self.ferw = np.zeros((nk, self.nband, self.nspin_band))
patt = re.compile("^orbital") patt = re.compile("^orbital")
# FIXME: fix spin indices for NCDIJ = 4 (non-collinear)
assert self.ncdij < 4, "Non-collinear case is not implemented"
for ispin in range(self.nspin): for ispin in range(self.nspin):
for ik in range(nk): for ik in range(nk):
for ib in range(self.nband): for ib in range(self.nband):
@ -301,7 +241,7 @@ class Plocar:
ctmp = complex(float(sline[1]), float(sline[2])) ctmp = complex(float(sline[1]), float(sline[2]))
plo[ip, ispin, ik, ib] = ctmp plo[ip, ispin, ik, ib] = ctmp
print("Read parameters:") print("Read parameters: LOCPROJ")
for il, par in enumerate(proj_params): for il, par in enumerate(proj_params):
print(il, " -> ", par) print(il, " -> ", par)
@ -504,7 +444,7 @@ class Kpoints:
self.kpts[ik, :] = list(map(float, sline[:3])) self.kpts[ik, :] = list(map(float, sline[:3]))
self.kwghts[ik] = float(sline[3]) self.kwghts[ik] = float(sline[3])
self.kwghts /= np.sum(self.kwghts) self.kwghts /= self.nktot
# Attempt to read tetrahedra # Attempt to read tetrahedra
# Skip comment line ("Tetrahedra") # Skip comment line ("Tetrahedra")
@ -636,7 +576,6 @@ class Doscar:
# First line: NION, NION, JOBPAR, NCDIJ # First line: NION, NION, JOBPAR, NCDIJ
sline = next(f).split() sline = next(f).split()
self.ncdij = int(sline[3])
# Skip next 4 lines # Skip next 4 lines
for _ in range(4): for _ in range(4):

13
python/triqs_dft_tools/converters/vasp.py
View File

@ -186,12 +186,11 @@ class VaspConverter(ConverterTools):
raise "VaspConverter: error reading %s"%self.ctrl_file raise "VaspConverter: error reading %s"%self.ctrl_file
# if nc_flag: # if nc_flag:
## TODO: check this # VASP.6.
# n_spin_blocs = 1 if SO == 1:
# else: n_spin_blocs = 1
# n_spin_blocs = ns else:
n_spin_blocs = SP + 1 - SO n_spin_blocs = SP + 1
# Read PLO groups # Read PLO groups
# First, we read everything into a temporary data structure # First, we read everything into a temporary data structure
# TODO: think about multiple shell groups and how to map them on h5 structures # TODO: think about multiple shell groups and how to map them on h5 structures
@ -283,7 +282,7 @@ class VaspConverter(ConverterTools):
# raise NotImplementedError("Noncollinear calculations are not implemented") # raise NotImplementedError("Noncollinear calculations are not implemented")
# else: # else:
hopping = numpy.zeros([n_k, n_spin_blocs, nb_max, nb_max], numpy.complex_) hopping = numpy.zeros([n_k, n_spin_blocs, nb_max, nb_max], numpy.complex_)
f_weights = numpy.zeros([n_k, n_spin_blocs, nb_max], numpy.float_) f_weights = numpy.zeros([n_k, n_spin_blocs, nb_max], numpy.complex_)
band_window = [numpy.zeros((n_k, 2), dtype=int) for isp in range(n_spin_blocs)] band_window = [numpy.zeros((n_k, 2), dtype=int) for isp in range(n_spin_blocs)]
n_orbitals = numpy.zeros([n_k, n_spin_blocs], numpy.int) n_orbitals = numpy.zeros([n_k, n_spin_blocs], numpy.int)