diff --git a/python/triqs_dft_tools/converters/plovasp/elstruct.py b/python/triqs_dft_tools/converters/plovasp/elstruct.py index 6418b5ea..cafa4d45 100644 --- a/python/triqs_dft_tools/converters/plovasp/elstruct.py +++ b/python/triqs_dft_tools/converters/plovasp/elstruct.py @@ -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,41 +150,75 @@ 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) + 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] - 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] - 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 -# ov_max = np.maximum(ov, ov_max) -# ov_min = np.minimum(ov, ov_min) + 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) diff --git a/python/triqs_dft_tools/converters/plovasp/plotools.py b/python/triqs_dft_tools/converters/plovasp/plotools.py index 1dba65de..50e064d9 100644 --- a/python/triqs_dft_tools/converters/plovasp/plotools.py +++ b/python/triqs_dft_tools/converters/plovasp/plotools.py @@ -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]) diff --git a/python/triqs_dft_tools/converters/plovasp/proj_group.py b/python/triqs_dft_tools/converters/plovasp/proj_group.py index 0237ae6e..71bb8601 100644 --- a/python/triqs_dft_tools/converters/plovasp/proj_group.py +++ b/python/triqs_dft_tools/converters/plovasp/proj_group.py @@ -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] @@ -417,7 +417,7 @@ class ProjectorGroup: # Calculate [O^{-1/2}]_{m m'} eig, eigv = np.linalg.eigh(overlap) assert np.all(eig > 0.0), ("Negative eigenvalues of the overlap matrix:" - "projectors are ill-defined") + "projectors are ill-defined") sqrt_eig = 1.0 / np.sqrt(eig) shalf = np.dot(eigv * sqrt_eig, eigv.conj().T) # Apply \tilde{P}_{m v} = \sum_{m'} [O^{-1/2}]_{m m'} P_{m' v} diff --git a/python/triqs_dft_tools/converters/plovasp/proj_shell.py b/python/triqs_dft_tools/converters/plovasp/proj_shell.py index f334f264..32461f57 100644 --- a/python/triqs_dft_tools/converters/plovasp/proj_shell.py +++ b/python/triqs_dft_tools/converters/plovasp/proj_shell.py @@ -115,7 +115,11 @@ class ProjectorShell: to avoid superfluous matrix multiplications. """ 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: 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 - nlm = self.lm2 - self.lm1 - _, ns, nk, nb = proj_raw.shape + nproj, 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: -# 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,14 +267,22 @@ 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(isp, ns_band) + 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[:, 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 +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,11 +373,13 @@ class ProjectorShell: el_struct.eigvals[:, ib1:ib2, isp]): for io in range(nion): proj_k = self.proj_win[io, isp, ik, ...] - 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) + # 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 return loc_ham diff --git a/python/triqs_dft_tools/converters/plovasp/vaspio.py b/python/triqs_dft_tools/converters/plovasp/vaspio.py index 97e1d6f6..9ccc67e8 100644 --- a/python/triqs_dft_tools/converters/plovasp/vaspio.py +++ b/python/triqs_dft_tools/converters/plovasp/vaspio.py @@ -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: @@ -256,6 +182,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") @@ -271,19 +205,25 @@ class Plocar: proj_params[ip]['label'] = label proj_params[ip]['isite'] = isite 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() - + assert ip == nproj, "Number of projectors in the header is wrong in LOCPROJ" self.eigs = np.zeros((nk, self.nband, self.nspin_band)) 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): @@ -299,9 +239,9 @@ 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:") + 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,8 +576,7 @@ 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): sline = next(f) diff --git a/python/triqs_dft_tools/converters/vasp.py b/python/triqs_dft_tools/converters/vasp.py index 5b5fa478..93f770a7 100644 --- a/python/triqs_dft_tools/converters/vasp.py +++ b/python/triqs_dft_tools/converters/vasp.py @@ -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)