Merge ../bare_plovasp into dev

Conflicts:
	python/converters/vasp/python/main.py
	python/converters/vasp/python/plotools.py
Merged old version of output with the new one

python/converters/vasp/doc/source/plotools.rst

@@ -72,7 +72,7 @@ Consistency with parameters
 
 
 Selecting projector subsets
----------------------------
+===========================
 
 The first step of PLO processing is to select subsets of projectors
 corresponding to PLO groups. Each group contains a set of shells.
@@ -117,7 +117,7 @@ The class is using a helper function `select_bands()` for selecting a subset of
 .. _ortho:
 
 Orthogonalization
-=================
+-----------------
 
 At the second stage the selected projectors are orthogonalized (orthonormalized). 
 Orthogonalization can be performed in different ways if projection is made
@@ -138,3 +138,54 @@ The way the normalization of a PLO group is done is controlled by two group para
     (False by default)
 
 
+Storing generated projectors
+****************************
+
+After the PLOs are generated they are stored to text files which can be subsequently
+converted to TRIQS h5-files (using the converter). The general format of the file
+is a JSON-header containing all necessary parameters followed by a set of arrays.
+There is always one (control) file containing general information (`k`-kpoints, lattice vectors etc.)
+and `at least` one file containing correlated groups (one file for each group).
+
+Control file format
+===================
+
+Filename '<namebase>.ctrl'. Contains the data shared between all shells.
+The JSON-header consists of the following elements:
+
+  * *nk*: number of `k`-points
+
+  * *ns*: number of spin channels
+
+  * *nc_flag*: collinear/noncollinear case (False/True)
+
+  * *ng*: number of projector groups
+
+  * Symmetry information (list of symmetry operations)
+
+  * *efermi*: Fermi level (optional)
+
+  * Lattice information
+
+Projector-group file format
+===========================
+
+Projector-group files have names '<namebase>.plog<Ng>'.
+They essentially contain serialized objects of class 'ProjectorGroup'.
+The JSON-header has, thus, the following elements:
+
+  * *shells*: list of shells
+
+  * each shell is a dictionary:
+
+    - *lshell*: orbital number `l`
+
+    - *nion*: number of ions
+
+    - *ndim*: number of orbitals/ion
+
+    - *nbmax*: maxmimum number of bands (needed for array allocations)
+
+  * *emin*, *emax*: energy window
+
+

python/converters/vasp/python/inpconf.py

@@ -68,6 +68,9 @@ class ConfigParameters:
             'normalize' : ('normalize', self.parse_string_logical, True),
             'normion' : ('normion', self.parse_string_logical, False)}
 
+        self.gen_optional = {
+            'basename' : ('basename', str, 'vasp'),
+            'efermi' : ('efermi', float)}
 
 #
 # Special parsers
@@ -177,7 +180,7 @@ class ConfigParameters:
             key = param_set[par][0]
             try:
                 par_str = self.cp.get(section, par)
-            except ConfigParser.NoOptionError:
+            except (ConfigParser.NoOptionError, ConfigParser.NoSectionError):
                 if exception:
                     message = "Required parameter '%s' not found in section [%s]"%(par, section)
                     raise Exception(message)
@@ -438,8 +441,18 @@ class ConfigParameters:
         """
         Parses [General] section.
         """
-# TODO: write the parser
-        pass
+        self.general = {}
+        sections = self.cp.sections()
+        gen_section = filter(lambda s: s.lower() == 'general', sections)
+# If no [General] section is found parse a dummy section name to the parser
+# to reset parameters to their default values
+        if len(gen_section) > 1:
+            raise Exception("More than one section [General] is found")
+        if len(gen_section) == 0:
+            gen_section = 'general'
+        gen_section = gen_section[0]
+        parsed = self.parse_parameter_set(gen_section, self.gen_optional, exception=False)
+        self.general.update(parsed)
 
 ################################################################################
 #

python/converters/vasp/python/main.py

@@ -3,7 +3,7 @@ import sys
 import vaspio
 from inpconf import ConfigParameters
 from elstruct import ElectronicStructure
-from plotools import generate_plo, plo_output, kpoints_output
+from plotools import generate_ortho_plos, output_as_text
 
 if __name__ == '__main__':
     narg = len(sys.argv)
@@ -23,11 +23,5 @@ if __name__ == '__main__':
     pars.parse_input()
     vasp_data = vaspio.VaspData(vasp_dir)
     el_struct = ElectronicStructure(vasp_data)
-    pshells, pgroups = generate_plo(pars, el_struct)
-    for gr in pgroups:
-        gr.orthogonalize()
-
-# TODO: add BASENAME to config parameters
-    basename = 'vasp'
-    kpoints_output(basename, el_struct)
-    plo_output(basename, pshells, pgroups, el_struct)
+    pshells, pgroups = generate_ortho_plos(pars, el_struct)
+    output_as_text(pars, el_struct, pshells, pgroups)

python/converters/vasp/python/plotools.py

@@ -2,6 +2,9 @@
 import itertools as it
 import numpy as np
 
+# 'simplejson' is supposed to be faster than 'json' in stdlib.
+import simplejson as json
+
 class Projector:
     """
     Class describing a local-orbital projector.
@@ -277,6 +280,12 @@ class ProjectorShell:
         self.lm2 = (self.lorb+1)**2
         self.ndim = self.lm2 - self.lm1
 
+        if self.tmatrix is None:
+            self.ndim = self.lm2 - self.lm1 + 1
+        else:
+# TODO: generalize this to a tmatrix for every ion
+            self.ndim = self.tmatrix.shape[0]
+
 # Pre-select a subset of projectors (this should be an array view => no memory is wasted)
 # !!! This sucks but I have to change the order of 'ib' and 'ilm' indices here
 # This should perhaps be done right after the projector array is read from PLOCAR
@@ -424,12 +433,49 @@ def kpoints_output(basename, el_struct):
 
 
 
+################################################################################
+#
+# ctrl_output
+#
+################################################################################
+def ctrl_output(conf_pars, el_struct, ng):
+    """
+    Outputs a ctrl-file.
+    """
+    ctrl_fname = conf_pars.general['basename'] + '.ctrl'
+    head_dict = {}
+
+# TODO: Add output of tetrahedra
+# Construct the header dictionary
+    head_dict['ngroups'] = ng
+    head_dict['nk'] = el_struct.kmesh['nktot']
+    head_dict['ns'] = el_struct.nspin
+    head_dict['nc_flag'] = 1 if el_struct.nc_flag else 0
+#    head_dict['efermi'] = conf_pars.general['efermi']  # We probably don't need Efermi
+
+    header = json.dumps(head_dict, indent=4, separators=(',', ': '))
+
+    print "  Storing ctrl-file..."
+    with open(ctrl_fname, 'wt') as f:
+        f.write(header + "\n")
+        f.write("#END OF HEADER\n")
+
+        f.write("# k-points and weights\n")
+        labels = ['kx', 'ky', 'kz', 'kweight']
+        out = "".join(map(lambda s: s.center(15), labels)) 
+        f.write("#" + out + "\n")
+        for ik, kp in enumerate(el_struct.kmesh['kpoints']):
+            tmp1 = "".join(map("{0:15.10f}".format, kp))
+            out = tmp1 + "{0:16.10f}".format(el_struct.kmesh['kweights'][ik])
+            f.write(out + "\n")
+
+
 ################################################################################
 #
 # plo_output
 #
 ################################################################################
-def plo_output(basename, pshells, pgroups, el_struct):
+def plo_output(conf_pars, el_struct, pshells, pgroups):
     """
     Outputs PLO groups into text files.
 
@@ -455,55 +501,91 @@ def plo_output(basename, pshells, pgroups, el_struct):
     # Projected shells
     Nshells
     # Shells: <shell indices>
+    # Shell <1>
     Shell 1
     ndim
     # complex arrays: plo(ns, nion, ndim, nb)
     ...
     # Shells: <shell indices>
+    # Shell <2>
     Shell 2
     ...
 
     """
-    for ig, gr in enumerate(pgroups):
-        fname = basename + '.plog%i'%(ig+1)
-        with open(fname, 'wt') as f:
-            f.write("# Energy window: emin, emax\n")
-            f.write("%i   %i\n"%(gr.emin, gr.emax))
-            f.write("# Number of electrons within the window\n")
-            nelect = gr.nelect_window(el_struct)
-            f.write("%s\n"%(round(nelect, 5)))
-            f.write("# Eigenvalues: is, ik, ib1, ib2 then list of values\n")
+    for ig, pgroup in enumerate(pgroups):
+        plo_fname = conf_pars.general['basename'] + '.pg%i'%(ig + 1)
+        print "  Storing PLO-group file '%s'..."%(plo_fname)
+        head_dict = {}
 
-            nk, ns_band, _ = gr.ib_win.shape
+        head_dict['ewindow'] = (pgroup.emin, pgroup.emax)
+        head_dict['nb_max'] = pgroup.nb_max
+
+        head_shells = []
+        for ish in pgroup.ishells:
+            shell = pgroup.shells[ish]
+            sh_dict = {}
+            sh_dict['shell_index'] = ish
+            sh_dict['lorb'] = shell.lorb
+            sh_dict['ndim'] = shell.ndim
+# Convert ion indices from the internal representation (starting from 0)
+# to conventional VASP representation (starting from 1)
+            ion_output = [io + 1 for io in shell.ion_list]
+            sh_dict['ion_list'] = ion_output
+# TODO: add the output of transformation matrices
+
+            head_shells.append(sh_dict)
+
+        head_dict['shells'] = head_shells
+
+        header = json.dumps(head_dict, indent=4, separators=(',', ': '))
+
+        with open(plo_fname, 'wt') as f:
+            f.write(header + "\n")
+            f.write("#END OF HEADER\n")
+            
+# Eigenvalues within the window
+            f.write("# Eigenvalues within the energy window: %s, %s\n"%(pgroup.emin, pgroup.emax))
+            nk, nband, ns_band = el_struct.eigvals.shape
             for isp in xrange(ns_band):
+                f.write("# is = %i\n"%(isp + 1))
                 for ik in xrange(nk):
-                    ib1 = gr.ib_win[ik, isp, 0]
-                    ib2 = gr.ib_win[ik, isp, 1]
-                    f.write("%i  %i    %i  %i\n"%(isp+1, ik+1, ib1+1, ib2+1))
+                    ib1, ib2 = pgroup.ib_win[ik, isp, 0], pgroup.ib_win[ik, isp, 1]
+                    f.write(" %i  %i\n"%(ib1, ib2))
                     for ib in xrange(ib1, ib2 + 1):
-                        f.write("%s\n"%(el_struct.eigvals[ik, ib, isp]))
+                        f.write("%15.7f\n"%(el_struct.eigvals[ik, ib, isp]))
 
-            f.write("\nProjected shells: nshells\n")
-            f.write("%i\n"%(len(gr.ishells)))
-            f.write("Shells: <shell indices>\n")
-            f.write('  '.join(map(lambda ish: "{0:d}".format(pshells[ish].user_index), gr.ishells)) + '\n')
-            for ish in gr.ishells:
-                shell = pshells[ish]
-                f.write("Orbital dimension: ndim\n")
-                f.write("%i\n"%(shell.ndim))
+# Projected shells
+            f.write("# Projected shells\n")
+            f.write("# Shells: %s\n"%(pgroup.ishells))
+            for ish in pgroup.ishells:
+                shell = pgroup.shells[ish]
+                f.write("# Shell %i\n"%(ish))
 
-                nion, ns, nk, ndim, _ = shell.proj_win.shape
-                f.write("# Blocks [isp, ion, ndim, nb], 'nb' fastest index\n")
-                for ik in xrange(nk):
+                nion, ns, nk, nlm, nb = shell.proj_win.shape
                 for isp in xrange(ns):
-# TODO: fix this for non-collinear case (ns != ns_band)
-                        ib1 = gr.ib_win[ik, isp, 0] - gr.nb_min
-                        ib2 = gr.ib_win[ik, isp, 1] - gr.nb_min
+                    f.write("# is = %i\n"%(isp + 1))
+                    for ik in xrange(nk):
+                        f.write("# ik = %i\n"%(ik + 1))
                         for ion in xrange(nion):
-                            f.write("# ik = %i, is = %i, ion = %i\n"%(ik+1, isp+1, ion+1))
-                            for iorb in xrange(ndim):
-                                for ib in xrange(ib1, ib2+1):
-                                    plo = shell.proj_win[ion, isp, ik, iorb, ib]
-                                    f.write("%20.10f %20.10f\n"%(plo.real, plo.imag))
+                            for ilm in xrange(nlm):
+                                ib1, ib2 = pgroup.ib_win[ik, isp, 0], pgroup.ib_win[ik, isp, 1]
+                                ib1_win = ib1 - shell.nb_min
+                                ib2_win = ib2 - shell.nb_min
+                                for ib in xrange(ib1_win, ib2_win + 1):
+                                    p = shell.proj_win[ion, isp, ik, ilm, ib]
+                                    f.write("{0:16.10f}{1:16.10f}\n".format(p.real, p.imag))
                                 f.write("\n")
 
+
+################################################################################
+#
+# output_as_text
+#
+################################################################################
+def output_as_text(pars, el_struct, pshells, pgroups):
+    """
+    Output all information necessary for the converter as text files.
+    """
+    ctrl_output(pars, el_struct, len(pgroups))
+    plo_output(pars, el_struct, pshells, pgroups)
+

python/converters/vasp/test/_inpconf/.gitignore

@@ -0,0 +1 @@
+*.pyc

python/converters/vasp/test/_inpconf/example.cfg

@@ -1,4 +1,6 @@
 [General]
+BASENAME = test_base
+EFERMI = 0.1
 
 [Group 1]
 SHELLS = 1 2

python/converters/vasp/test/_inpconf/test_general.py

@@ -0,0 +1,32 @@
+r"""
+Tests of 'parse_general()' defined in ConfigParameters class
+"""
+import arraytest
+import numpy as np
+from inpconf import ConfigParameters
+
+################################################################################
+#
+# TestParseGeneral
+#
+################################################################################
+class TestParseGeneral(arraytest.ArrayTestCase):
+    """
+    Function:
+
+    def parse_general(self)
+
+    Scenarios:
+
+    - **if** a correct [General] section is defined **return** a dictionary
+    """
+# Scenario 1
+    def test_example(self):
+        conf_pars = ConfigParameters('example.cfg')
+        conf_pars.parse_general()
+        res = conf_pars.general
+        expected = {'basename': 'test_base', 'efermi': 0.1}
+        self.assertDictEqual(res, expected)
+
+
+

python/converters/vasp/test/adv_example.cfg

@@ -1,20 +0,0 @@
-[General]
-EFERMI = -0.6
-
-[Group 1]
-SHELLS = 1 2
-EMIN = -7.6
-EMAX =  2.7
-
-[Shell 1]
-# Ni shell
-IONS = 5 6 7 8
-LSHELL = 2
-RTRANSFORM =
-  0.0 0.0 0.0 0.0 1.0
-  0.0 0.0 1.0 0.0 0.0
-
-[Shell 2]
-# Oxygen shell
-IONS = 9..20
-LSHELL = 1

python/converters/vasp/test/run_all.sh

@@ -0,0 +1 @@
+PYTHONPATH=$HOME/Codes/vasp/vasp5.3/plo_vasp/plovasp/python:$HOME/Codes/vasp/vasp5.3/plo_vasp/plovasp/c:$PYTHONPATH python test_all.py

python/converters/vasp/test/test_all.py

@@ -0,0 +1,10 @@
+r"""
+Searches and runs all available test suites.
+"""
+import unittest
+
+if __name__ == '__main__':
+    suite = unittest.TestLoader().discover('./')
+    unittest.TextTestRunner(verbosity=2, buffer=True).run(suite)
+#    unittest.TextTestRunner(verbosity=2, buffer=False).run(suite)
+