WIP: Vasp RPA optics

README.rst

@@ -117,9 +117,12 @@ Feature support for different codes
 +------------------+----------+--------+----------+
 | phonon band plot | **Y** (2)| **Y**  |   *N*    |
 +------------------+----------+--------+----------+
+| optical spectra  | **Y** (3)|  *N*   |  **Y**   |
++------------------+----------+--------+----------+
 
 (1) Brillouin-zone path can also be written for CASTEP phonon calculation
 (2) VASP phonons are plotted from Phonopy output files
+(3) As well as the LOPTICS sum over empty states, RPA spectra with ALGO=CHI are supported
 
 Installation
 ------------

docs/source/sumo-optplot.rst

@@ -1,9 +1,9 @@
 sumo-optplot
 ==============
 
-``sumo-optplot`` is a program for generating publication-ready optical absorption
-spectra diagrams from VASP calculations. The script supports plotting multiple
-spectra simultaneously.
+``sumo-optplot`` is a program for generating publication-ready optical
+absorption spectra diagrams from VASP or Questaal calculations. The
+script supports plotting multiple spectra simultaneously.
 
 .. contents:: Table of Contents
    :local:
@@ -18,7 +18,7 @@ and be can be accessed using the command::
     sumo-optplot -h
 
 To plot an absorption spectra, simply run the following command in a folder containing a ``vasprun.xml`` or
-``vasprun.xml.gz`` file, which has been calculated using ``LOPTICS = .TRUE.``::
+``vasprun.xml.gz`` file, which has been calculated using ``LOPTICS = .TRUE.`` or ``ALGO = CHI``::
 
     sumo-optplot
 

sumo/cli/optplot.py

@@ -19,7 +19,7 @@
 
 from pymatgen.io.vasp import Vasprun
 from pymatgen.util.string import latexify
-from sumo.io import questaal
+from sumo.io import questaal, vasp_rpa
 
 from sumo.plotting.optics_plotter import SOpticsPlotter
 from sumo.electronic_structure.optics import (broaden_eps,
@@ -105,7 +105,7 @@ def optplot(modes=('absorption',), filenames=None, codes='vasp',
 
     ##### BUILD LIST OF FILES AUTOMATICALLY IF NECESSARY #####
 
-    if codes == 'vasp':
+    if codes in ('vasp', 'vasp-rpa'):
         if not filenames:
             if os.path.exists('vasprun.xml'):
                 filenames = ['vasprun.xml']
@@ -153,6 +153,11 @@ def optplot(modes=('absorption',), filenames=None, codes='vasp',
             else:
                 auto_band_gaps.append(None)
 
+        elif code == 'vasp-rpa':
+            dielectrics.append(
+                vasp_rpa.dielectric_from_vasprun(filename))
+            auto_band_gaps.append(None)
+
         elif code == 'questaal':
             if not save_files:
                 out_filename = None

sumo/io/vasp_rpa.py

@@ -0,0 +1,69 @@
+"""Read dielectric data from Vasp RPA calculations"""
+
+from xml.etree import ElementTree
+import numpy as np
+
+def dielectric_from_vasprun(filename, theory='auto'):
+    """Read a vasprun.xml file and return dielectric function from RPA
+
+    This function can interpret calculations performed with ALGO=CHI. For
+    regular LOPTICS Vasp runs, use the standard (Pymatgen-based) importer
+
+    Args:
+        filename (:obj:`float`):
+            Path to ``vasprun.xml`` output of ALGO=CHI calculation.
+        theory (:obj:`str`):
+            Calculation type to use. Depending on LRPA tag, 'rpa' or 'lfe'
+            (local field effects including XC field) is available. If 'auto',
+            use whichever of these methods is available. In either case 'ipa'
+            (independent particle approximation) is available.
+
+    Returns:
+        :obj:`tuple`
+            The format imitates the ``dielectric`` attribute of
+            :obj:`pymatgen.io.vasp.outputs.Vasprun`: a tuple of the form::
+
+                ([energy1, energy2, ...],
+                 [[real_xx_1, real_yy_1, real_zz_1,
+                   real_xy_1, real_yz_1, real_xz_1],
+                  [real_xx_2, real_yy_2, real_zz_2,
+                   real_xy_2, real_yz_2, real_xz_2], ...],
+                 [[imag_xx_1, imag_yy_1, imag_zz_1,
+                   imag_xy_1, imag_yz_1, imag_xz_1],
+                  [imag_xx_2, imag_yy_2, imag_zz_2,
+                   imag_xy_2, imag_yz_2, imag_xz_2], ...])
+    """
+
+    vrxml = ElementTree.parse(filename)
+    root = vrxml.getroot()
+    dielectric_functions = root.findall('dielectricfunction')    
+
+    calcs_tags = {'ipa': 'INDEPENDENT PARTICLE',
+                  'rpa': 'RPA',
+                  'lfe': 'local field effects'}
+
+    data = {}
+    for calc, tag in calcs_tags.items():
+        try:
+            diel_xml = next((diel for diel in dielectric_functions
+                            if tag in diel.attrib['comment']))
+            data[calc] = {}
+            for component in ('real', 'imag'):
+                text_rows = ((el.text for el in
+                              diel_xml.find(component).find('array').find('set')))
+                data[calc][component] = np.genfromtxt(text_rows)
+        except StopIteration:
+            continue
+
+    theory = theory.lower()
+    if theory == 'auto':
+        if 'rpa' in data:
+            theory = 'rpa'
+        elif 'lfe' in data:
+            theory = 'lfe'
+        else:
+            raise Exception("No RPA or LFE data in vasprun")
+
+    return (data[theory]['real'][:, 0],
+            data[theory]['real'][:, 1:],
+            data[theory]['imag'][:, 1:])