Exciting: added xas functionality

Ocean: Calculation of eps_11 eps_22 and eps_33 in one run is now possible

exciting_handler.py

@@ -11,6 +11,7 @@
 import re
 import threading
 from little_helpers import convert_to_ordered
+from sortedcontainers import SortedSet
 
 p_table = {i: el.__repr__() for i, el in enumerate(pt.elements)}
 p_table_rev = {el.__repr__(): i for i, el in enumerate(pt.elements)}
@@ -168,7 +169,7 @@ def __init__(self):
              'screentype': 'full'
                 , 'nempty_screening': '0', 'use gw': 'false', 'nempty': '5', 'ngridq': "4 4 4", 'ngridk': "4 4 4",
              'broad': '0.0005',
-             'gqmax': "3.0", 'vkloff': "0.231 0.151 0.432"})
+             'gqmax': "3.0", 'vkloff': "0.231 0.151 0.432", 'xas': 'false', 'xasspecies': '1','xasatom': '1', 'nstlxas': '1 6','xasedge':'K'})
         self.optical_spectrum_options_tooltip = {
             'nstlbse': 'Range of bands included for the BSE calculation.\nThe first pair of numbers corresponds to the band index for local orbitals and valence states (counted from the lowest eigenenergy),\nthe second pair corresponds to the band index of the conduction states (counted from the Fermi level).',
             'nempty_screening': 'Number of empty states.',
@@ -199,7 +200,25 @@ def __init__(self):
 IP (Independent particles)
 RPA (random phase approximation)
 singlet
-triplet"""}
+triplet""",
+            'xas':'Set to "true" to perform BSE X-rasy absorption spectroscopy (XAS) calculation',
+            'xasatom':'Atom number (as defined in the structure in OpenDFT) for which the XAS is calculated.',
+            'xasspecies':'Species number (enumerate the atoms in your structure from low Z to high Z) for which the XAS is calculated. Example LiBH_4: H:1 Li:2 B:3',
+            'xasedge':"""Defines the initial states of the XAS calculation.
+Type: 	choose from:
+K
+L1
+L2
+L3
+L23
+M1
+M2
+M3
+M23
+M4
+M5
+M45""",
+            'nstlxas':'Range of bands included for the BSE calculation. The pair corresponds to the band index of the conduction states (counted from the Fermi level).'}
 
         self.relax_file_timestamp = None
 
@@ -670,7 +689,7 @@ def _add_scf_to_tree(self, tree, crystal_structure, skip=False):
             ET.SubElement(crystal, "basevect").text = "{0:1.6f} {1:1.6f} {2:1.6f}".format(*lattice_vector)
 
         abs_coord_atoms = crystal_structure.atoms
-        species = set(abs_coord_atoms[:, 3].astype(np.int))
+        species = SortedSet(abs_coord_atoms[:, 3].astype(np.int))
         n_species = len(species)
         n_atoms = abs_coord_atoms.shape[0]
 
@@ -729,7 +748,9 @@ def _add_optical_spectrum_to_tree(self, tree):
         screening = ET.SubElement(xs, "screening", screentype=self.optical_spectrum_options['screentype'],
                                   nempty=self.optical_spectrum_options['nempty_screening'])
         BSE = ET.SubElement(xs, "BSE", bsetype=self.optical_spectrum_options['bsetype'],
-                            nstlbse=self.optical_spectrum_options['nstlbse'])
+                            nstlbse=self.optical_spectrum_options['nstlbse'],xas=self.optical_spectrum_options['xas'],
+                            xasspecies=self.optical_spectrum_options['xasspecies'],xasatom=self.optical_spectrum_options['xasatom'],
+                            xasedge=self.optical_spectrum_options['xasedge'],nstlxas=self.optical_spectrum_options['nstlxas'])
         qpointset = ET.SubElement(xs, 'qpointset')
         qpoint = ET.SubElement(qpointset, 'qpoint').text = '0.0 0.0 0.0'
 

ocean_handlers.py

@@ -7,7 +7,7 @@
 import solid_state_tools as sst
 import numpy as np
 import shutil
-
+from sortedcontainers import SortedSet
 
 
 class OceanHandler(object):
@@ -17,9 +17,10 @@ def __init__(self):
         self.working_dirctory = '/abinit_ocean_files/'
 
         self.optical_spectrum_options = convert_to_ordered({'diemac':'5.0','CNBSE.xmesh':'6 6 6','screen.shells':'4.0','screen.shells':'3.5','cnbse.rad':'3.5',
-        'cnbse.broaden':'0.1','edges':'0 1 0','screen.nbands':'80','screen.nkpt':'2 2 2','opts.core_states':'1 0 0 0','opts.relativity':'scalar rel','opts.functional':'lda','opts.valence':'2.0 3.5 0.0 0.0',
+        'cnbse.broaden':'0.1','edges':'0 1 0','screen.nbands':'80','screen.nkpt':'2 2 2','opts.core_states':'1 0 0 0',
+                                                            'opts.relativity':'scalar rel','opts.functional':'lda','opts.valence':'2.0 3.5 0.0 0.0',
         'fill.pow':'2','fill.energy':'0.30 2.00 0.0001','fill.cutoff':'3.5','fill.fourier':'0.05 20',
-        'operator':'dipole','polarization':'0 0 1','k vector':'0 1 0','energy range':'600'})
+        'operator':'dipole','polarization':'all','k vector':'0 1 0','energy range':'600'})
 
     def start_optical_spectrum(self, crystal_structure):
         """This method starts a optical spectrum calculation in a subprocess. The configuration is stored in optical_spectrum_options.
@@ -54,22 +55,34 @@ def read_optical_spectrum(self):
     - optical_spectrum:       A OpticalSpectrum object with the latest optical spectrum result found.
         """
 
-        files = []
-        for file in os.listdir(self.project_directory+self.working_dirctory+'/CNBSE'):
-            if file.startswith("absspct_") and file.endswith('_01'):
-                files.append(file)
+        def load_eps(i):
+            files = []
+            for file in os.listdir(self.project_directory+self.working_dirctory+'/CNBSE'):
+                if file.startswith("absspct_") and file.endswith('_0{}'.format(i+1)):
+                    files.append(file)
+
+            if not files:
+                return None,None,None
 
 
-        data = np.loadtxt(self.project_directory+self.working_dirctory+'/CNBSE/'+files[0], skiprows=2)
+            data = np.loadtxt(self.project_directory+self.working_dirctory+'/CNBSE/'+files[0], skiprows=2)
 
-        for file in files[1:]:
-            data += np.loadtxt(self.project_directory+self.working_dirctory+'/CNBSE/'+file, skiprows=2)
+            for file in files[1:]:
+                data += np.loadtxt(self.project_directory+self.working_dirctory+'/CNBSE/'+file, skiprows=2)
 
-        energy = data[:, 0]
-        eps2 = data[:, 1]
-        eps1 = data[:, 3]
+            energy = data[:, 0]
+            eps2 = data[:, 1]
+            eps1 = data[:, 3]
+            return energy,eps1,eps2
 
-        return sst.OpticalSpectrum(energy, eps2, epsilon1=eps1)
+        eps1_list = []
+        eps2_list = []
+        for i in range(3):
+            energy,eps1,eps2 = load_eps(i)
+            eps1_list.append(eps1)
+            eps2_list.append(eps2)
+
+        return sst.OpticalSpectrum(energy, eps2_list, epsilon1=eps1_list)
 
     def _make_ocean_input_file(self, filename='ocean.in'):
         if not os.path.isdir(self.project_directory + self.working_dirctory):
@@ -85,9 +98,9 @@ def _add_ocean_to_file(self,file,crystal_structure):
 
         file.write('nbands '+self.scf_options['nband']+'\n')
 
-        atom_species = set(crystal_structure.atoms[:,3])
+        atom_species = SortedSet(crystal_structure.atoms[:,3])
         file.write('pp_list{\n')
-        for specie in sorted(atom_species):
+        for specie in atom_species:
             file.write(p_table[specie]+'.fhi\n' )
         file.write('}\n')
 
@@ -141,17 +154,24 @@ def _make_fill_file(self,filename='ocean.fill'):
 
         f.close()
 
-    def _make_photon_file(self,filename='photon1'):
+    def _make_photon_file(self):
         if not os.path.isdir(self.project_directory + self.working_dirctory):
             os.mkdir(self.project_directory + self.working_dirctory)
-        f = open(self.project_directory + self.working_dirctory + '/' + filename, 'w')
 
-        f.write(self.optical_spectrum_options['operator']+'\n')
-        f.write('cartesian '+self.optical_spectrum_options['polarization']+'\n')
-        f.write('end\n')
-        f.write('cartesian '+self.optical_spectrum_options['k vector']+'\n')
-        f.write('end\n')
-        f.write(self.optical_spectrum_options['energy range'])
+        if self.optical_spectrum_options['polarization'].strip().lower() == 'all':
+            polarizations = ['1 0 0','0 1 0','0 0 1']
+        else:
+            polarizations = [self.optical_spectrum_options['polarization']]
+
+        for i,polarization in enumerate(polarizations):
+            f = open(self.project_directory + self.working_dirctory + '/' + 'photon{0}'.format(i+1), 'w')
+
+            f.write(self.optical_spectrum_options['operator']+'\n')
+            f.write('cartesian '+polarization+'\n')
+            f.write('end\n')
+            f.write('cartesian '+self.optical_spectrum_options['k vector']+'\n')
+            f.write('end\n')
+            f.write(self.optical_spectrum_options['energy range'])
 
         f.close()