From 3317371762357d60757a9a26161e430ff3aad22a Mon Sep 17 00:00:00 2001 From: "Oleg E. Peil" Date: Thu, 19 Feb 2015 21:22:51 +0100 Subject: [PATCH] Added calculation of density matrices for a shell A method 'density_matrix()' for evaluating a density matrix of a given shell has been added to class ProjectorShell. It requires an ElectronicStructure object as an input an by default produces a site- and spin-diagonal part of the density matrix using the Fermi-weights obtained directly from VASP. Ideally, this density matrix should coincide with the one calculated within VASP itself (inside the LDA+U module). Corresponding sanity test has been added, which shows only that the calculation does not crash. Real numerical tests are needed. --- python/converters/vasp/python/elstruct.py | 10 ++++-- python/converters/vasp/python/plotools.py | 34 +++++++++++++++++++ python/converters/vasp/python/vaspio.py | 14 ++++---- .../vasp/test/plotools/test_projgroups.py | 8 +++-- .../vasp/test/plotools/test_projshells.py | 30 ++++++++++++++++ 5 files changed, 86 insertions(+), 10 deletions(-) diff --git a/python/converters/vasp/python/elstruct.py b/python/converters/vasp/python/elstruct.py index dc378e69..7a9f3cc0 100644 --- a/python/converters/vasp/python/elstruct.py +++ b/python/converters/vasp/python/elstruct.py @@ -15,6 +15,7 @@ class ElectronicStructure: - *efermi* (float) : Fermi level read from DOSCAR - *proj_raw* (array[complex]) : raw projectors from PLOCAR - *eigvals* (array[float]) : KS eigenvalues + - *ferw* (array[float]) : Fermi weights from VASP - *kmesh* (dict) : parameters of the `k`-mesh - *structure* (dict) : parameters of the crystal structure - *symmetry* (dict) : paramters of symmetry @@ -29,12 +30,12 @@ class ElectronicStructure: self.kmesh = {'nktot': self.nktot} self.kmesh['kpoints'] = vasp_data.kpoints.kpts - self.kmesh['kwights'] = vasp_data.eigenval.kwghts + self.kmesh['kweights'] = vasp_data.eigenval.kwghts try: self.kmesh['ntet'] = vasp_data.kpoints.ntet self.kmesh['itet'] = vasp_data.kpoints.itet self.kmesh['volt'] = vasp_data.kpoints.volt - except NameError: + except AttributeError: pass # Note that one should not subtract this Fermi level from eigenvalues @@ -55,6 +56,11 @@ class ElectronicStructure: # [see ProjectorGroup.orthogonalization()] self.proj_raw = vasp_data.plocar.plo.transpose((0, 1, 2, 4, 3)) +# In fact, Fermi weights do not depend on ions +# FIXME: restructure the data in PLOCAR to remove the redundant dependency +# of 'ferw' on ions + self.ferw = vasp_data.plocar.ferw[0, :, :, :] + # Check that the number of atoms is the same in PLOCAR and POSCAR natom_plo = vasp_data.plocar.params['nion'] assert natom_plo == self.natom, "PLOCAR is inconsistent with POSCAR (number of atoms)" diff --git a/python/converters/vasp/python/plotools.py b/python/converters/vasp/python/plotools.py index f46311b0..5655b9ea 100644 --- a/python/converters/vasp/python/plotools.py +++ b/python/converters/vasp/python/plotools.py @@ -1,4 +1,5 @@ +import itertools as it import numpy as np class Projector: @@ -292,6 +293,39 @@ class ProjectorShell: ib2_win = ib2 - self.nb_min self.proj_win[:, isp, ik, :, ib1_win:ib2_win] = self.proj_arr[:, isp, ik, :, ib1:ib2] +################################################################################ +# +# select_projectors +# +################################################################################ + def density_matrix(self, el_struct, site_diag=True, spin_diag=True): + """ + Returns occupation matrix/matrices for the shell. + """ + nion, ns, nk, nlm, nbtot = self.proj_win.shape + + assert site_diag, "site_diag = False is not implemented" + assert spin_diag, "spin_diag = False is not implemented" + + occ_mats = np.zeros((ns, nion, nlm, nlm), dtype=np.float64) + + kweights = el_struct.kmesh['kweights'] + occnums = el_struct.ferw + ib1 = self.nb_min + ib2 = self.nb_max + 1 + for isp in xrange(ns): + for ik, weight, occ in it.izip(it.count(), kweights, occnums[isp, :, :]): + for io in xrange(nion): + proj_k = self.proj_win[isp, io, ik, ...] + occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2], + 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 + + def generate_ortho_plos(conf_pars, vasp_data): """ diff --git a/python/converters/vasp/python/vaspio.py b/python/converters/vasp/python/vaspio.py index 2a22464a..f3d96105 100644 --- a/python/converters/vasp/python/vaspio.py +++ b/python/converters/vasp/python/vaspio.py @@ -54,7 +54,9 @@ class Plocar: Properties ---------- - plo (numpy.array((nion, ns, nk, nb, nlmmax))) : raw projectors + - *plo* (numpy.array((nion, ns, nk, nb, nlmmax))) : raw projectors + - *params* (dict) : parameters read from PLOCAR + - *ferw* (array(nion, ns, nk, nb)) : Fermi weights from VASP """ def from_file(self, vasp_dir='./', plocar_filename='PLOCAR'): @@ -206,11 +208,11 @@ class Kpoints: Properties ---------- - nktot (int) : total 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 + - nktot (int) : total 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) """ # diff --git a/python/converters/vasp/test/plotools/test_projgroups.py b/python/converters/vasp/test/plotools/test_projgroups.py index 26e0c607..de238537 100644 --- a/python/converters/vasp/test/plotools/test_projgroups.py +++ b/python/converters/vasp/test/plotools/test_projgroups.py @@ -46,7 +46,9 @@ class TestProjectorGroup(mytest.MyTestCase): ib1 = self.proj_gr.ib_win[ik, 0, 0] ib2 = self.proj_gr.ib_win[ik, 0, 1] f.write("%i %i\n"%(ib1, ib2)) - for ib in xrange(ib2 - self.proj_gr.nb_min + 1): + ib1w = ib1 - self.proj_gr.nb_min + ib2w = ib2 - self.proj_gr.nb_min + 1 + for ib in xrange(ib1w, ib2w): for ilm in xrange(nlm): p = self.proj_gr.shells[0].proj_win[ion, isp, ik, ilm, ib] f.write("%5i %s\n"%(ilm+1, p)) @@ -65,7 +67,9 @@ class TestProjectorGroup(mytest.MyTestCase): ib1 = self.proj_gr.ib_win[ik, 0, 0] ib2 = self.proj_gr.ib_win[ik, 0, 1] f.write("%i %i\n"%(ib1, ib2)) - for ib in xrange(ib2 - self.proj_gr.nb_min + 1): + ib1w = ib1 - self.proj_gr.nb_min + ib2w = ib2 - self.proj_gr.nb_min + 1 + for ib in xrange(ib1w, ib2w): for ilm in xrange(nlm): p = self.proj_gr.shells[0].proj_win[ion, isp, ik, ilm, ib] f.write("%5i %s\n"%(ilm+1, p)) diff --git a/python/converters/vasp/test/plotools/test_projshells.py b/python/converters/vasp/test/plotools/test_projshells.py index ebb81799..c78bd7f5 100644 --- a/python/converters/vasp/test/plotools/test_projshells.py +++ b/python/converters/vasp/test/plotools/test_projshells.py @@ -1,6 +1,7 @@ import numpy as np import vaspio +import elstruct from inpconf import ConfigParameters from plotools import select_bands, ProjectorShell import mytest @@ -18,6 +19,7 @@ class TestProjectorShell(mytest.MyTestCase): Scenarios: - compare output for a correct input + - test density matrix """ # Scenario 1 def test_example(self): @@ -53,4 +55,32 @@ class TestProjectorShell(mytest.MyTestCase): expected_file = 'projshells.out' self.assertFileEqual(testout, expected_file) + +# Scenario 2 + def test_dens_mat(self): + conf_file = 'example.cfg' + pars = ConfigParameters(conf_file) + pars.parse_input() + vasp_data = vaspio.VaspData('./') + el_struct = elstruct.ElectronicStructure(vasp_data) + + efermi = el_struct.efermi + eigvals = el_struct.eigvals - efermi + emin = pars.groups[0]['emin'] + emax = pars.groups[0]['emax'] + ib_win, nb_min, nb_max = select_bands(eigvals, emin, emax) + + proj_sh = ProjectorShell(pars.shells[0], vasp_data.plocar.plo) + + proj_sh.select_projectors(ib_win, nb_min, nb_max) + + dens_mat = proj_sh.density_matrix(el_struct) + print dens_mat + + testout = 'densmat.out.test' + with open(testout, 'wt') as f: + f.write("density matrix: %s\n"%(dens_mat)) + + expected_file = 'densmat.out' + self.assertFileEqual(testout, expected_file)