mirror of
https://github.com/triqs/dft_tools
synced 2024-12-22 20:34:38 +01:00
Merge pull request #218 from dariofiosca/unstable
Non-collinear projectors VASP.6.
This commit is contained in:
Alexander Hampel
GitHub
commit
aab83893cf
@ -61,6 +61,8 @@ class ElectronicStructure:
|
||||
|
||||
self.kmesh = {'nktot': self.nktot}
|
||||
self.kmesh['kpoints'] = vasp_data.kpoints.kpts
|
||||
# VASP.6.
|
||||
self.nc_flag = vasp_data.plocar.nc_flag
|
||||
self.kmesh['kweights'] = vasp_data.kpoints.kwghts
|
||||
try:
|
||||
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
|
||||
# of bands in case of non-collinear calculations
|
||||
self.nspin = vasp_data.plocar.nspin
|
||||
self.nc_flag = vasp_data.plocar.ncdij == 4
|
||||
|
||||
self.nband = vasp_data.plocar.nband
|
||||
|
||||
# Check that the number of k-points is the same in all files
|
||||
@ -150,12 +150,11 @@ class ElectronicStructure:
|
||||
norb = nproj // nions
|
||||
|
||||
# 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
|
||||
|
||||
if self.nc_flag == False:
|
||||
den_mat = np.zeros((ns, nproj, nproj), dtype=np.float64)
|
||||
overlap = np.zeros((ns, nproj, nproj), dtype=np.float64)
|
||||
# ov_min = np.ones((ns, nproj, nproj), dtype=np.float64) * 100.0
|
||||
# ov_max = np.zeros((ns, nproj, nproj), dtype=np.float64)
|
||||
for ispin in range(ns):
|
||||
for ik in range(nk):
|
||||
kweight = self.kmesh['kweights'][ik]
|
||||
@ -163,9 +162,6 @@ class ElectronicStructure:
|
||||
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
|
||||
# ov_max = np.maximum(ov, ov_max)
|
||||
# ov_min = np.minimum(ov, ov_min)
|
||||
|
||||
# Output only the site-diagonal parts of the matrices
|
||||
print()
|
||||
print(" Unorthonormalized density matrices and overlaps:")
|
||||
@ -188,3 +184,41 @@ 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=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)
|
||||
|
||||
|
||||
|
||||
@ -85,7 +85,7 @@ def check_data_consistency(pars, el_struct):
|
||||
# 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"
|
||||
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])
|
||||
|
||||
@ -111,7 +111,7 @@ class ProjectorGroup:
|
||||
"""
|
||||
self.nelect = 0
|
||||
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 ik in range(nk):
|
||||
ib1 = self.ib_win[ik, isp, 0]
|
||||
|
||||
@ -115,7 +115,11 @@ class ProjectorShell:
|
||||
to avoid superfluous matrix multiplications.
|
||||
"""
|
||||
nion = self.nion
|
||||
# VASP.6.
|
||||
if self.nc_flag == False:
|
||||
nm = self.lm2 - self.lm1
|
||||
else:
|
||||
nm = 2*(self.lm2 - self.lm1)
|
||||
|
||||
if 'tmatrices' in sh_pars:
|
||||
self.do_transform = True
|
||||
@ -133,21 +137,7 @@ class ProjectorShell:
|
||||
nr = nrow // nion
|
||||
nsize = ncol // nm
|
||||
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
|
||||
ns_dim = max(1, nsize // 2)
|
||||
|
||||
@ -196,7 +186,7 @@ class ProjectorShell:
|
||||
# If no transformation matrices are provided define a default one
|
||||
self.do_transform = False
|
||||
|
||||
ns_dim = 2 if self.nc_flag else 1
|
||||
ns_dim = 1
|
||||
ndim = nm * ns_dim
|
||||
|
||||
# We still need the matrices for the output
|
||||
@ -219,16 +209,15 @@ class ProjectorShell:
|
||||
according to the shell parameters.
|
||||
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
|
||||
nproj, ns, nk, nb = proj_raw.shape
|
||||
if self.nc_flag == 0:
|
||||
nlm = self.lm2 - self.lm1
|
||||
_, ns, nk, nb = proj_raw.shape
|
||||
else:
|
||||
nlm = 2*(self.lm2 - self.lm1)
|
||||
|
||||
if self.do_transform:
|
||||
# TODO: implement a non-collinear case
|
||||
# for a non-collinear case 'ndim' is 'ns * nm'
|
||||
ndim = self.tmatrices.shape[1]
|
||||
self.proj_arr = np.zeros((nion, ns, nk, ndim, nb), dtype=np.complex128)
|
||||
for ik in range(nk):
|
||||
@ -236,8 +225,6 @@ class ProjectorShell:
|
||||
for io, ion in enumerate(self.ion_list):
|
||||
proj_k = np.zeros((ns, nlm, nb), dtype=np.complex128)
|
||||
qcoord = structure['qcoords'][ion]
|
||||
# kphase = np.exp(-2.0j * np.pi * np.dot(kp, qcoord))
|
||||
# kphase = 1.0
|
||||
for m in range(nlm):
|
||||
# Here we search for the index of the projector with the given isite/l/m indices
|
||||
for ip, par in enumerate(proj_params):
|
||||
@ -257,14 +244,8 @@ class ProjectorShell:
|
||||
for ip, par in enumerate(proj_params):
|
||||
if par['isite'] - 1 == ion and par['l'] == self.lorb and par['m'] == m:
|
||||
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
|
||||
|
||||
|
||||
################################################################################
|
||||
#
|
||||
# select_projectors
|
||||
@ -286,9 +267,17 @@ class ProjectorShell:
|
||||
|
||||
# Select projectors for a given energy window
|
||||
ns_band = self.ib_win.shape[1]
|
||||
for isp in range(ns):
|
||||
if ns == 1:
|
||||
for ik in range(nk):
|
||||
# TODO: for non-collinear case something else should be done here
|
||||
is_b = min(0, 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[:, :, 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
|
||||
@ -322,15 +311,26 @@ class ProjectorShell:
|
||||
occnums = el_struct.ferw
|
||||
ib1 = self.ib_min
|
||||
ib2 = self.ib_max + 1
|
||||
# VASP.6.
|
||||
count_nan = 0
|
||||
print("Site diag : {}".format(site_diag))
|
||||
if site_diag:
|
||||
for isp in range(ns):
|
||||
for ik, weight, occ in zip(it.count(), kweights, occnums[isp, :, :]):
|
||||
for io in range(nion):
|
||||
proj_k = self.proj_win[io, isp, ik, ...]
|
||||
occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2],
|
||||
proj_k.conj().T).real * weight
|
||||
overlaps[isp, io, :, :] += np.dot(proj_k,
|
||||
proj_k.conj().T).real * weight
|
||||
# VASP.6.
|
||||
array_sum = np.sum(proj_k)
|
||||
if np.isnan(array_sum) == True:
|
||||
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:
|
||||
proj_k = np.zeros((ndim, nbtot), dtype=np.complex128)
|
||||
for isp in range(ns):
|
||||
@ -339,13 +339,11 @@ class ProjectorShell:
|
||||
i1 = io * nlm
|
||||
i2 = (io + 1) * nlm
|
||||
proj_k[i1:i2, :] = self.proj_win[io, isp, ik, ...]
|
||||
# TODO
|
||||
occ_mats[isp, 0, :, :] += np.dot(proj_k * occ[ib1:ib2],
|
||||
proj_k.conj().T).real * weight
|
||||
overlaps[isp, 0, :, :] += np.dot(proj_k,
|
||||
proj_k.conj().T).real * weight
|
||||
overlaps[isp, 0, :, :] += np.dot(proj_k,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
|
||||
|
||||
@ -375,12 +373,14 @@ class ProjectorShell:
|
||||
el_struct.eigvals[:, ib1:ib2, isp]):
|
||||
for io in range(nion):
|
||||
proj_k = self.proj_win[io, isp, ik, ...]
|
||||
# VASP.6.
|
||||
if self.nc_flag == True:
|
||||
loc_ham[isp, io, :, :] += np.dot(proj_k * (eigk - el_struct.efermi),
|
||||
proj_k.conj().T) * weight * 0.5
|
||||
else:
|
||||
loc_ham[isp, io, :, :] += np.dot(proj_k * (eigk - el_struct.efermi),
|
||||
proj_k.conj().T) * weight
|
||||
|
||||
# if not symops is None:
|
||||
# occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map)
|
||||
|
||||
return loc_ham
|
||||
|
||||
################################################################################
|
||||
|
||||
@ -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.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'):
|
||||
r"""
|
||||
@ -242,8 +164,12 @@ class Plocar:
|
||||
line = f.readline()
|
||||
line = line.split("#")[0]
|
||||
sline = line.split()
|
||||
self.ncdij, nk, self.nband, nproj = list(map(int, sline[:4]))
|
||||
self.nspin = 1 if self.ncdij == 1 else 2
|
||||
self.ncdij, nk, self.nband, nproj = list(map(int, sline[0:4]))
|
||||
|
||||
# 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
|
||||
|
||||
try:
|
||||
@ -257,6 +183,14 @@ class Plocar:
|
||||
iproj_site = 0
|
||||
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
|
||||
line = self.search_for(f, "^ *ISITE")
|
||||
ip = 0
|
||||
@ -271,9 +205,16 @@ class Plocar:
|
||||
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
|
||||
else:
|
||||
proj_params[ip]['m'] = 2*m + 1
|
||||
else:
|
||||
proj_params[ip]['m'] = m
|
||||
|
||||
ip +=1
|
||||
|
||||
line = f.readline().strip()
|
||||
|
||||
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))
|
||||
|
||||
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 ik in range(nk):
|
||||
for ib in range(self.nband):
|
||||
@ -301,7 +241,7 @@ class Plocar:
|
||||
ctmp = complex(float(sline[1]), float(sline[2]))
|
||||
plo[ip, ispin, ik, ib] = ctmp
|
||||
|
||||
print("Read parameters:")
|
||||
print("Read parameters: LOCPROJ")
|
||||
for il, par in enumerate(proj_params):
|
||||
print(il, " -> ", par)
|
||||
|
||||
@ -504,7 +444,7 @@ class Kpoints:
|
||||
self.kpts[ik, :] = list(map(float, sline[:3]))
|
||||
self.kwghts[ik] = float(sline[3])
|
||||
|
||||
self.kwghts /= np.sum(self.kwghts)
|
||||
self.kwghts /= self.nktot
|
||||
|
||||
# Attempt to read tetrahedra
|
||||
# Skip comment line ("Tetrahedra")
|
||||
@ -636,7 +576,6 @@ class Doscar:
|
||||
|
||||
# First line: NION, NION, JOBPAR, NCDIJ
|
||||
sline = next(f).split()
|
||||
self.ncdij = int(sline[3])
|
||||
|
||||
# Skip next 4 lines
|
||||
for _ in range(4):
|
||||
|
||||
@ -186,12 +186,11 @@ class VaspConverter(ConverterTools):
|
||||
raise "VaspConverter: error reading %s"%self.ctrl_file
|
||||
|
||||
# if nc_flag:
|
||||
## TODO: check this
|
||||
# n_spin_blocs = 1
|
||||
# else:
|
||||
# n_spin_blocs = ns
|
||||
n_spin_blocs = SP + 1 - SO
|
||||
|
||||
# VASP.6.
|
||||
if SO == 1:
|
||||
n_spin_blocs = 1
|
||||
else:
|
||||
n_spin_blocs = SP + 1
|
||||
# Read PLO groups
|
||||
# First, we read everything into a temporary data structure
|
||||
# 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")
|
||||
# else:
|
||||
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)]
|
||||
n_orbitals = numpy.zeros([n_k, n_spin_blocs], numpy.int)
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user