def display(string):
  '''Prints :literal:`string` to the terminal output and to the .oklog file.'''
  global log_fid
  if not options.quiet:
    print(string)
  if not options.no_log:
    if log_fid is None:
      if options.outputname is None or options.outputname == '':
        return
      else:
        log_fid = '%s.oklog' % options.outputname.split('@')[0]
    if not is_initiated:
      init_display(log_fid)
    f = open(log_fid, 'a')
    f.write('%s\n' % string)
    f.close() class QCinfo:
  '''Class managing all information from the from the output
  files of quantum chemical software.

  See :ref:`Central Variables` in the manual for details.
  '''
  def __init__(self, data=None):
    self.geo_info = []
    self.geo_spec = []
    
    if data:
      self.geo_spec = data['geo_spec']
      self.geo_info = data['geo_info']
      self.format_geo()
      if isinstance(data['ao_spec'], numpy.ndarray):
        ao_spec = data['ao_spec'][numpy.newaxis][0]
      else:
        ao_spec = data['ao_spec']
      if isinstance(data['mo_spec'], numpy.ndarray):
        mo_spec = data['mo_spec'][numpy.newaxis][0]
      else:
        mo_spec = data['mo_spec']
      self.ao_spec = AOClass(restart=ao_spec)
      self.mo_spec = MOClass(restart=mo_spec)
    else:
      self.geo_info = []
      self.geo_spec = []
      self.ao_spec = AOClass()
      self.mo_spec = MOClass()

  def __eq__(self, other):
    if not isinstance(other, QCinfo):
      raise TypeError('Comaring of QCinfo to non QCinfo object not defined')
    same = [
            self.comp_geo_info(other.geo_info),
            numpy.allclose(self.geo_spec, other.geo_spec),
            self.ao_spec == other.ao_spec,
            self.mo_spec == other.mo_spec
            ]
    return all(same)

  def update(self):
    self.ao_spec.update()
    self.mo_spec.update()

  def comp_geo_info(self, geo2):
    same = True
    for atom1, atom2 in zip(self.geo_info, geo2):
      if not len(atom1) == len(atom2):
        raise ValueError('Atom object are of different length!')
      for i in range(len(atom1)):
        if atom1[i] != atom2[i]:
          same = False
    return same

  def copy(self):
    from copy import deepcopy
    qcinfo = deepcopy(self)
    qcinfo.update()
    return qcinfo

  def format_geo(self, is_angstrom=False):
    '''Converts geo_info and geo_spec to a universal format.
    **Parameters:**

    is_angstrom : bool, optional
      If True, input is assumed to be in Angstrom and positions are converted to Bohr radii.
    '''
    for i in self.geo_info:
      i[0] = get_atom_symbol(i[0])
      i[2] = float(i[-1])
    self.geo_info = numpy.array(self.geo_info)
    self.geo_spec = numpy.array(self.geo_spec,dtype=float)
    if is_angstrom:
      self.geo_spec *= aa_to_a0

  def get_com(self,nuc_list=None):
    '''Computes the center of mass.
    '''
    com   = numpy.zeros(3)
    total_mass = 0.
    if nuc_list is None:
      nuc_list = list(range(len(self.geo_spec))) # iterate over all nuclei
    for ii in nuc_list:
      nuc_mass    = standard_mass(self.geo_info[ii][0])
      com   += numpy.multiply(self.geo_spec[ii],nuc_mass)
      total_mass += nuc_mass
    com = com/total_mass
    return com


  def get_charge(self, nuclear=True, electron=True):
    '''Computes total charge of the system.
    '''
    charge = 0.
    if electron:
      charge -= sum(self.mo_spec.get_occ())
    if nuclear:
      for ii in range(len(self.geo_info)):
        charge += float(self.geo_info[ii][2])
    return charge

  def get_coc(self):
    '''Computes the center of charge.
    '''
    coc     = numpy.zeros(3)
    for ii in range(len(self.geo_info)):
      nuc_charge    = float(self.geo_info[ii][2])
      coc += numpy.multiply(self.geo_spec[ii],nuc_charge)
    coc = coc / self.get_charge(nuclear=True)
    return coc

  def get_bc(self,matrix=None,is_vector=False):
    '''Calculates Barycenter for scalar field
    '''
    # Initialize variable
    self.bc = numpy.zeros(3)
    # Calculation of barycenter
    from orbkit import grid
    if not is_vector:
      grid.grid2vector()
    xyz = grid.tolist()
    for i in range(3):
      self.bc[i] = (matrix.reshape((-1,))*xyz[i]).sum()
    self.bc /= matrix.sum()
    if not is_vector:
      grid.vector2grid(*grid.N_)
    return self.bc

  def select_spin(self,restricted,spin=None):
    '''For an unrestricted calculation, the name of the MO
    ('sym' keyword in ``qc.mo_spec``) is modified, e.g., 3.1_b for MO 3.1 with
    beta spin and 3.1_a for MO 3.1 with alpha spin.
    For restricted calculation, the 'spin' keyword from ``qc.mo_spec`` is
    removed.

    **Parameters:**

    restricted : bool
      If True, removes the 'spin' keyword from ``qc.mo_spec``.
    spin : {None, 'alpha', or 'beta'}, optional
      If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
    '''
    # Only molecular orbitals of one spin requested?
    if spin is not None:
      for i in range(len(self.mo_spec))[::-1]:
        if self.mo_spec[i]['spin'] != spin:
          del self.mo_spec[i]

    if restricted:
      # Closed shell calculation
      for mo in self.mo_spec:
        del mo['spin']
    else:
      # Rename MOs according to spin
      for mo in self.mo_spec:
        mo['sym'] += '_%s' % mo['spin'][0]
      if not isinstance(self.mo_spec, MOClass):
        self.mo_spec = MOClass(self.mo_spec)
        self.mo_spec.get_spinstate()

  def todict(self):
    '''Returns the dictionary that is used to save QCinfo instance
    '''
    data = {}
    data['ao_spec'] = self.ao_spec.todict()
    data['mo_spec'] = self.mo_spec.todict()
    data['geo_spec'] = self.geo_spec
    data['geo_info'] = self.geo_info
    data['parent_class_name'] = self.__module__ + '.' + self.__class__.__name__
    return data

  def get_ase_atoms(self,bbox=None,**kwargs):
    '''Create an ASE atoms object.
    (cf. https://wiki.fysik.dtu.dk/ase/ase/atoms.html )

    **Parameters:**

    bbox : list of floats (bbox=[xmin,xmax,ymin,ymax,zmin,zmax]), optional
      If not None, sets the unit cell to the grid boundaries and moves the
      molecule in its center.

    **Returns:**

    atoms : Atoms object
      See https://wiki.fysik.dtu.dk/ase/ase/atoms.html for details

    .. Note::

      ASE has to be in the PYTHONPATH
    '''
    from ase import Atoms
    from ase.units import Bohr

    atoms = Atoms("".join(self.geo_info[:,0]),
                  positions=self.geo_spec*Bohr,
                  **kwargs)
    if bbox is not None:
      if len(bbox) != 6:
        raise ValueError("bbox has to have 6 elements")
      bbox = numpy.array(bbox)
      atoms.translate(-bbox[::2]*Bohr)
      atoms.cell = numpy.eye(3) * (bbox[1::2] - bbox[::2])*Bohr

    return atoms
  # Synonym
  atoms = get_ase_atoms

  def view(self,select=slice(None,None,None),bbox=None,**kwargs):
    '''Opens ase-gui with the atoms of the QCinfo class.
    (cf. https://wiki.fysik.dtu.dk/ase/ase/visualize/visualize.html )

    **Parameters:**

    select : slice or (array of int), default: all atoms
      Specifies the atoms to be shown.
    bbox : list of floats (bbox=[xmin,xmax,ymin,ymax,zmin,zmax]), optional
      If not None, sets the unit cell to the grid boundaries and moves the
      molecule in its center.

    .. Note::

      ASE has to be in the PYTHONPATH
    '''
    from ase import visualize
    visualize.view(self.get_ase_atoms(bbox=bbox,**kwargs)[select])

  @property
  def nuclear_repulsion(self):
    '''Calculates nuclear repulsion energy.'''
    from numpy.linalg import norm
    Vnn = 0
    Natoms = self.geo_spec.shape[0]
    for a in range(Natoms):
      Za = float(self.geo_info[a,2])
      Ra = self.geo_spec[a,:].astype(float)
      for b in range(a+1, Natoms):
        Zb = float(self.geo_info[b,2])
        Rb = self.geo_spec[b,:].astype(float)
        Vnn += Za*Zb / norm(Ra-Rb)
    return Vnn class CIinfo():
  '''Class managing all information from the from the output
  files of quantum chemical software for CI calculations.

  The CI related features are in ongoing development.
  '''
  def __init__(self,method='ci'):
    self.method = method
    self.info   = None
    self.coeffs = []
    self.occ    = []
    self.moocc  = None

  def __str__(self):
    string = '%s' % self.method.upper()
    if self.info is not None:
      string += ' State %(state)s' % self.info
      if 'spin' in self.info.keys() and self.info['spin'] != 'Unknown':
         string += ' (%(spin)s)' % self.info
    if numpy.shape(self.coeffs) != (0,):
      string += ':\tNorm = %0.8f (%d Coefficients)' %(self.get_norm(),
                                                      len(self.coeffs))
    return  string
  def __eq__(self, other):
    try:
      return self.__dict__ == other.__dict__
    except ValueError:
      return False
  def get_norm(self):
    return sum(self.coeffs**2)
  def renormalize(self):
    self.coeffs /= self.get_norm()
  def apply_threshold(self,threshold,keep_length=False):
    i = numpy.abs(self.coeffs) > threshold
    if keep_length:
      self.coeffs[numpy.invert(i)] = 0.0
    else:
      if self.info['state'] == '0':
       self.coeffs = self.coeffs[i]
       self.occ = self.occ[:1]
      else:
       self.coeffs = self.coeffs[i]
       self.occ = self.occ[i]
      
  def copy(self):
    ciinfo = deepcopy(self)
    return ciinfo

  def get_moocc(self):
    if self.moocc is None:
      raise ValueError('ci.set_moocc(qc) has to be called first! (ci.moocc is not initialized)')
    return self.moocc

  def set_moocc(self,moocc):
    assert (moocc.dtype == numpy.intc), 'moocc has to be numpy.intc'
    self.moocc = moocc

  def hdf5_save(self,fid='out.h5',group='/ci:0',mode='w'):
    from orbkit.output import hdf5_open,hdf5_append
    from copy import copy
    for hdf5_file in hdf5_open(fid,mode=mode):
      dct = copy(self.todict())
      dct['info'] = numpy.array(dct['info'].items(),dtype=str)
      hdf5_append(dct,hdf5_file,name=group)

  def hdf5_read(self,fid='out.h5',group='/ci:0'):
    from orbkit.output import hdf5_open,hdf52dict
    for hdf5_file in hdf5_open(fid,mode='r'):
      for key in self.__dict__.keys():
        try:
          self.__dict__[key] = hdf52dict('%s/%s' % (group,key),hdf5_file)
        except KeyError:
          self.__dict__[key] = hdf5_file['%s' % group].attrs[key]
      self.__dict__['info'] = dict(self.__dict__['info']) class MOClass(UserList):
  '''MO base class which contains all information on atomic orbitals.
  Two types of dataformats are available:

  1. Numpy-style data:

      coeffs : numpy.ndarray, dtype=float64, shape = (NMO, NAO)
        Molecular orbital coefficients.
      occ : numpy.ndarray, dtype=float64, shape = (NMO)
        Occupation numbers for molecular orbitals.
      eig : numpy.ndarray, dtype=float64, shape = (NMO)
        Eigenvalues for molecular orbitals.
      sym : numpy.ndarray, dtype=str, shape = (NMO)
        MOLPRO-like symmetry label of molecular orbitals.

  2. Lists of dictionaries:

    mo_spec

    See :ref:`Central Variables` in the manual for details.
  '''
  def __init__(self, seq = [], restart=None, data = None):
    if isinstance(data, list):
      seq = data
    elif isinstance(data, dict):
      restart = data
    UserList.__init__(self, seq)
    self._up_to_date = False
    self.coeffs = None
    self.occ = None
    self.sym = None
    self.spin = None
    self.eig = None
    self.spinpolarized = False
    self.selection_string = None
    self.selected_mo = None
    self.alpha_index = None
    self.beta_index = None
    if restart is not None:
      self._up_to_date = True
      self.coeffs = restart['coeffs']
      self.occ = restart['occ']
      self.sym = restart['sym']
      self.eig = restart['eig']
      self.spinpolarized = restart['spinpolarized']
      self.selection_string = restart.get('selection_string', None)
      self.selected_mo = restart.get('selected_mo', None)
      self.alpha_index = restart['alpha_index']
      self.beta_index = restart['beta_index']
      try:
        self.spin = restart['spin']
      except KeyError:
        warnings.warn('"spin" not found in restart file', UserWarning)
        self.spin = numpy.zeros_like(self.sym,dtype=str)
        self.spin[self.alpha_index] = 'alpha'
        self.spin[self.beta_index] = 'beta'
      self.internal_to_dict()

  def todict(self):
    self.update()
    data = {'coeffs': self.coeffs,
            'selection_string': self.selection_string,
            'alpha_index': self.alpha_index,
            'beta_index': self.beta_index,
            'selected_mo': self.selected_mo,
            'spinpolarized': self.spinpolarized,
            'occ': self.occ,
            'eig': self.eig,
            'parent_class_name' : self.__module__ + '.' + self.__class__.__name__,
            'sym': self.sym,
            'spin': self.spin}
    return data

  def __getitem__(self, item):
    parse_directly = False
    if isinstance(item, (int, numpy.int64)):
      return UserList.__getitem__(self, item)
    elif isinstance(item, slice):
      return MOClass(UserList.__getitem__(self, item))
    elif isinstance(item, (list, numpy.ndarray)) or \
         (sys.version_info.major == 3 and isinstance(item, range)):
      if isinstance(item, numpy.ndarray):
        parse_directly = item.dtype in [int, numpy.int_, numpy.intc, bool, numpy.bool_]
      else:
        parse_directly = True
        for it in item:
          if not isinstance(it, (int, numpy.int_, numpy.intc, bool, numpy.bool_)):
            parse_directly = False
            break
    if parse_directly:
      item = numpy.array(item)
      if item.ndim > 1:
        raise ValueError('Only 1D arrays can be used for indexing!')
      data_out = []
      for i, c in enumerate(item):
        if isinstance(c, numpy.bool_):
          if c:
            data_out.append(self.data[i])
        else:
          data_out.append(self.data[c])
      mo_out = MOClass(data_out)
      mo_out.selected_mo = item
      mo_out.update()
      del data_out
    else:
      mo_out = self.select(item)
    return mo_out

  def __setitem__(self, i, item):
    self.data[i] = item
    self._up_to_date = False

  def __delitem__(self, i):
    del self.data[i]
    self._up_to_date = False

  def __eq__(self, other):
    cases = [isinstance(other, MOClass), other == [], other is None]
    if not any(cases):
      raise TypeError('Comaring of MOClass to non MOClass object not defined')
    if cases[0]:
      self.update()
      same = [numpy.allclose(self.coeffs, other.coeffs),
      numpy.allclose(self.eig, other.eig),
      numpy.allclose(self.occ, other.occ),
      self.compare_sym(other.sym)]
      return all(same)
    else:
      if self.data is None or len(self.data) == 0:
        return True
      else:
        return False

  def compare_sym(self, sym2):
    same = True
    for atom1, atom2 in zip(self.sym, sym2):
      if not len(atom1) == len(atom2):
        raise ValueError('sym object are of different length!')
      for i in range(len(self.sym)):
        if self.sym[i] != sym2[i]:
          same = False
    return same

  def __str__(self):
    return '\n'.join(self.get_labels())

  def splinsplit_array(self, array):
    array_alpha = array[self.alpha_index]
    array_beta = array[self.beta_index]
    return array_alpha, array_beta

  @property
  def is_energy_sorted(self):
    return all(numpy.diff(self.get_eig())>=0)
  
  def get_homo(self, tol=1e-5, sort=True):
    '''Returns index of highest occupied MO.
    '''
    if not self._up_to_date:
      self.update()
    if sort:
      self.sort_by_energy()
    
    ihomo = (self.get_occ() > tol).nonzero()[0]
    
    return  None if not len(ihomo) else ihomo[-1]

  def get_lumo(self, tol=1e-5, sort=True):
    '''Returns index of lowest unoccupied MO.
    '''
    if not self._up_to_date:
      self.update()
    if sort:    
      self.sort_by_energy()
    
    ilumo = (self.get_occ() < tol).nonzero()[0]
    
    return None if not len(ilumo) else ilumo[0]

  def get_lastbound(self, sort=True):
    '''Returns index of highest bound MO.
    '''
    if not self._up_to_date:
      self.update()
    if sort:
      self.sort_by_energy()
    
    imaxbound = (self.get_eig() <= 0.).nonzero()[0]
    
    return None if not len(imaxbound) else imaxbound[-1]

  def sort_by_sym(self):
    '''Sorts mo_spec by symmetry.
    '''
    self.data.sort()
    self.update()

  def append(self, item):
    UserList.append(self, item)
    self._up_to_date = False

  def extend(self, item):
    UserList.extend(self, item)
    self._up_to_date = False

  def remove(self, item):
    UserList.remove(self, item)
    self._up_to_date = False

  def mo_template(self):
    template = {'coeffs': None,
                'energy': None,
                'occ_num': None,
                'sym': None,
                'spin': None}
    return template

  def internal_to_dict(self):
    self.data = []
    for imo in range(len(self.occ)):
      self.data.append(self.mo_template())
      self.data[-1]['coeffs'] = self.coeffs[imo]
      self.data[-1]['energy'] = self.eig[imo]
      self.data[-1]['occ_num'] = self.occ[imo]
      self.data[-1]['sym'] = self.sym[imo]
      self.data[-1]['spin'] = self.spin[imo]
    return

  def update(self):
    self._up_to_date = False
    self.get_coeffs()
    self.get_occ()
    self.get_eig()
    self.get_sym()
    self.get_spinstate()
    self._up_to_date = True
    return

  def sort_by_energy(self):
    if self.is_energy_sorted:
      return
    warnings.warn('MOs are not sorted by energy. Sorting them...', UserWarning)
    tmp_data = []
    sort = numpy.argsort(self.get_eig())
    for s in sort:
      tmp_data.append(self.data[s])
    self.data = tmp_data
    self.update()
    return

  def get_spinstate(self):
    '''Determines whether the MOClass has alpha and beta spins and removes the _a/_b spin labels.
    For spin-paired calculations all spins are set to 'alpha'.
    '''
    if not self._up_to_date:
      self.get_sym()
      
    self.spinpolarized = False
    self.alpha_index = []
    self.beta_index = []
    spindic = {'a': self.alpha_index, 
               'b': self.beta_index, 
               'u': self.alpha_index} # unknown is counted as alpha_index
    spinlabel = {'a':'alpha','b':'beta'}
    for isym,jsym in enumerate(self.sym):
      if jsym.endswith('_a') or jsym.endswith('_b'):
        spindic[jsym[-1]].append(isym)
        self.sym[isym] = jsym[:-2]
        self.spin[isym] = spinlabel[jsym[-1]]
      else:
        spindic[self.spin[isym][0]].append(isym)
    
    self.internal_to_dict()
    if len(self.beta_index) != 0:
      self.spinpolarized = True

  def get_labels(self,format='default'):
    if format=='short':
      return ['%(sym)s' %
                  i for i in self.data]
    elif format=='print':
      return ['%(sym)s (Occ = %(occ_num).2f, E = %(energy)+.4f E_h)' %
                  i for i in self.data]
      
    elif format in ['cb','cube','vmd']:
      return ['%(sym)s,Occ=%(occ_num).1f,E=%(energy)+.2f' %
                  i for i in self.data]
    else:
      return ['%(sym)s, Occ=%(occ_num).2f, E=%(energy)+.4f E_h' %
                    i for i in self.data]

  def get_indices(self):
    if self.selected_mo is not None:
      return self.selected_mo
    else:
      return list(range(len(self.data)))
  
  def set_coeffs(self, item):
    '''Set function for numpy array version of molecular orbital coefficients.

       **Returns:**

        coeff: numpy.ndarray, dtype=numpy.float64, shape = (NMO, NAO)
    '''
    require(item, dtype=numpy.float64)
    if not self.coeffs.shape == item.shape:
      raise ValueError('Old and new arrays need to be of the same size!')
    self.coeffs = item
    self.internal_to_dict()

  def set_occ(self, item):
    '''Set function for numpy array version of molecular orbital occupancies.

       **Parameters:**

        occ: numpy.ndarray, dtype=numpy.float64, shape = (NMO)
    '''
    require(item, dtype=numpy.float64)
    if not self.occ.shape == item.shape:
      raise ValueError('Old and new arrays need to be of the same size!')
    self.occ = item
    self.internal_to_dict()

  def set_eig(self, item):
    '''Set function for numpy array version of molecular orbital energies.

       **Parameters:**

        eig: numpy.ndarray, dtype=numpy.float64, shape = (NMO)
    '''
    require(item, dtype=numpy.float64)
    if not self.eig.shape == item.shape:
      raise ValueError('Old and new arrays need to be of the same size!')
    self.eig = item
    self.internal_to_dict()

  def set_sym(self, item):
    '''Set function for numpy array version of molecular orbital symmetries.

       **Parameters:**

        sym: numpy.ndarray, dtype=numpy.str, shape = (NMO)
    '''
    require(item, dtype=str)
    if not self.sym.shape == item.shape:
      raise ValueError('Old and new arrays need to be of the same size!')
    self.sym = item
    self.internal_to_dict()

  def get_coeffs(self):
    '''Get function for numpy array version of molecular orbital coefficients.

       **Returns:**

        coeffs: numpy.ndarray, dtype=numpy.float64, shape = (NMO, NAO)
    '''
    if not self._up_to_date:
      self.coeffs = numpy.zeros(shape=(len(self.data), len(self.data[0]['coeffs'])), dtype=numpy.float64)
      for imo, mo in enumerate(self.data):
        self.coeffs[imo] = mo['coeffs']
    return copy(self.coeffs)

  def get_eig(self):
    '''Get function for numpy array version of molecular orbital energies.

       **Returns:**

        eig: numpy.ndarray, dtype=numpy.float64, shape = (NMO)
    '''
    if not self._up_to_date:
      self.eig = numpy.zeros(shape=(len(self.data)), dtype=numpy.float64)
      for imo, mo in enumerate(self.data):
        self.eig[imo] = mo['energy']
    return copy(self.eig)

  def get_occ(self, return_alpha_beta=False, return_int=False, tol_int=1e-5, sum_occ=False):
    '''Get function for numpy array version of molecular orbital occupancies.

       **Parameters:**

        return_alpha_beta: bool, optional, determines whether occupancies should be
        returned as a shape = (NMO) array (False) or a shape = (2,NMO/2) array (True)
        return_int: bool, optional, determines whether occupancies should be
        returned as a numpy.intc array. Raises an error if there are partially occupied orbitals.
        tol_int: float, optional, defines the tolerance for when an occupation number is considered an integer
        sum_occ: bool, optional, determines whether occupancies should be summed over

       **Note:**

        return_alpha_beta=True is only supported for spin-polarized calculations.

       **Returns:**

        occ: numpy.ndarray, dtype=numpy.float64, shape = (NMO)
    '''
    if not self._up_to_date:
      self.occ = numpy.zeros(shape=(len(self.data)), dtype=numpy.float64)
      for imo, mo in enumerate(self.data):
        self.occ[imo] = mo['occ_num']
    if return_int:
      for f in self.occ:
        if self.spinpolarized and 0. > f < 1.-tol_int:
          raise ValueError('Occupation numbers are not integers')
        elif not self.spinpolarized and 0. > f < 2.-tol_int:
          raise ValueError('Occupation numbers are not integers')
    if sum_occ:
      return sum(self.get_occ(return_int=True))
    else:
      if not return_alpha_beta or not self.spinpolarized:
        if not return_int:
          return copy(self.occ)
        else:
          return copy(numpy.array(self.occ, dtype=numpy.intc))
      else:
        occ_alpha = self.occ[self.alpha_index]
        occ_beta = self.occ[self.beta_index]
        if not return_int:
          return copy(numpy.array([occ_alpha,occ_beta]))
        else:
          occ_alpha = numpy.array(occ_alpha, dtype=numpy.intc)
          occ_beta = numpy.array(occ_beta, dtype=numpy.intc)
          return copy(numpy.array([occ_alpha,occ_beta]))

  def get_sym(self):
    '''Get function for numpy array version of molecular orbital symmetries.

       **Returns:**

        sym: numpy.ndarray, dtype=numpy.str, shape = (NMO)
    '''
    if not self._up_to_date:
      self.sym = []
      self.spin = []
      for imo, mo in enumerate(self.data):
        self.sym.append(mo['sym'])
        if 'spin' in mo.keys():
          self.spin.append(mo['spin'])
        else:
          self.spin.append('unknown')
      
      self.sym = numpy.array(self.sym, dtype=str)
      
    return copy(self.sym)

  def get_spin_index(self, spin):
    '''Function used to select MO's by spin. A numpy.ndarray is returned
       which contains the indexes of MO's of the selected spin.

      **Parameters:**

      spin : 'str', can be either 'alpha' or 'beta'

      **Returns:**

      alpha : numpy.ndarray, dtype=numpy.intc

    '''
    spindic = {'alpha': self.alpha_index, 'beta': self.beta_index}
    return numpy.array(spindic[spin], dtype=numpy.intc)

  def combine_restricted_orbitals(self):
    '''Combines unrestricted orbitals with same MO coefficients to restricted orbitals.
    (For example Psi4 creates a molden file with unrestricted orbitals after a ROHF calculation.)
    '''

    if not self.spinpolarized:
      # orbitals are already restricted
      return

    # assert same energy
    ener_a, ener_b = self.eig[self.alpha_index], self.eig[self.beta_index]
    assert numpy.all(ener_a == ener_b), "MO eigenvalues do not match"

    # assert same MO coefficients
    coeff_a, coeff_b = self.coeffs[self.alpha_index], self.coeffs[self.beta_index]
    assert numpy.all(coeff_a == coeff_b), "MO coeffs do not match"

    # add occupation and select alpha orbitals only
    occ = self.occ[self.alpha_index] + self.occ[self.beta_index]
    self.data = self[self.alpha_index].data
    self.spinpolarized = False
    self.update()
    self.set_occ(occ)
    self.update()

  def select(self, fid_mo_list, flatten_input=True, sort_indices=True):
    '''Selects molecular orbitals from an external file or a list of molecular
       orbital labels.

    **Parameters:**

      mo_spec :
        See :ref:`Central Variables` for details.
      fid_mo_list : str, `'all_mo'`, or list
        | If fid_mo_list is a str, specifies the filename of the molecular orbitals list.
        | If fid_mo_list is 'all_mo', creates a list containing all molecular orbitals.
        | If fid_mo_list is a list, provides a list (or a list of lists) of molecular
          orbital labels.

      flatten_input : boolean, optional
        Specifies wheter lists of lists should be flattened so a single MOClass instance can be returned rather than a list of MOClass instances.

      sort_indices : boolean, optional
        Specifies wheter list of indexes should be sorted before it is returned. This is only supported if flatten_input is set to ``True``.


    **Supported Formats:**

      Integer List (Counting from **Zero**!)::

        1       2       3
        5       4
        homo    lumo+2:lumo+4

      List with Symmetry Labels::

        1.1     2.1     1.3
        1.1     4.1
        4.1     2.3     2.1

      ``alpha`` and ``beta`` can be used together with symmetry labels to restrict the selection to orbitals of that symmetry.
      This option is not supported for integer lists. Note also that ``alpha`` and ``beta`` only restrict selection within one
      string of inputs. If you which to spin-restrict orbitlas given as a list of strings please use ``all_alpha`` or ``all_beta``.
      You may also select all MO's of a given symmetry by specifying only the symmetry label
      with a preceeding asterisc, e.g. ``*.A1`` to get all orbitals of A1 symmetry. 

    **Returns:**

      List of MOClass instances containing the selected orbitals as well as further information on the selection criteria used
      If a sinlge list is used as in input and/or flatten_input=True, an MOClass instance is returned instead.
    '''
    import re
    display('\nProcessing molecular orbital list...')
    if flatten_input and isinstance(list(fid_mo_list)[0], (list, numpy.ndarray)):
      display('\nWarning! Flattening of input lists requested!')
      display('\nIf this was not intended please set ``flatten_input=False``')
    
    mo_in_file = []
    selected_mo = []
    sym_select = False

    def ordered_set(inlist):
      outlist = []
      for i in inlist:
        if i not in outlist:
          outlist.append(i)
      return outlist

    def eval_mp(i, r):
      if r not in i:
        return i
      else:
        tmp = i.split(r)
        if r == '-':
          i = int(tmp[0]) - int(tmp[1])
        elif r == '+':
          i = int(tmp[0]) + int(tmp[1])
        else:
          raise ArithmeticError('Unknown Operation in input string.')
        return str(i)

    def remove_empty(items):
      out = []
      for i, item in enumerate(items):
        if item and not item == 'None':
          out.append(item)
      return out

    def range_separators(items):
      pos_range = []
      for i, item in enumerate(items):
        if item == '#':
          pos_range.append(i)
      pos_range = numpy.array(pos_range)
      return pos_range

    def parse_nosym(item):
      keys = {'homo': str(self.get_homo()),
              'lumo': str(self.get_lumo()),
              'last_bound': str(self.get_lastbound()),
              'lastbound': str(self.get_lastbound())}
            
      warnings_str = {'homo': 'occupied',
                  'lumo': 'unoccupied',
                  'last_bound': 'bound'}
      
      for key in keys:
        if key in item and key in warnings_str.keys() and keys[key] == 'None':
          print('No {0} orbitals are present. Ignoring respective entries...'.format(warnings_str[key]))
        item = item.replace(key, keys[key])

      #This seems quite insane to me... I don't really know what else to do though...
      pos_range = range_separators(remove_empty(item.replace(':',':#:').split(':')))
      item = remove_empty(item.split(':'))
      if isinstance(item, list):
        for i in range(len(item)):
          for operation in ['-','+']:
            item[i] = eval_mp(item[i], operation)
      else:
        for operation in ['-','+']:
          item = eval_mp(item, operation)
      
      if len(pos_range) > 0:
        s = 1
        if numpy.allclose(pos_range,numpy.array([0])) and len(item) == 1: #Stuff like :b
          a = 0
          b = int(item[0])
        elif numpy.allclose(pos_range,numpy.array([1])) and len(item) == 1: #Stuff like a:
          a = int(item[0])
          b = len(self.data)
        elif numpy.allclose(pos_range,numpy.array([1])) and len(item) == 2: #Stuff like a:b
          a = int(item[0])
          b = int(item[1])
        elif numpy.allclose(pos_range,numpy.array([1,3])) and len(item) == 3: #Stuff like a:b:s
          a = int(item[0])
          b = int(item[1])
          s = int(item[2])
        elif numpy.allclose(pos_range,numpy.array([0,1])) and len(item) == 1: #Stuff like ::s
          a = 0
          b = len(self.data)
          s = int(item[0])
        else:
          raise ValueError('Format not recognized')
        if a < 0:
          a = 0
          warnings.warn('Lower index or index range out of bounds. Setting to zero', UserWarning)
        if b > len(self.data):
          b = len(self.data) - 1
          warnings.warn('Upper index or index range out of bounds. Setting to maximum index', UserWarning)
        item = [k for k in range(a, b, s)]
      else:
        item = [int(k) for k in item]
      return item

    def parse_sym(item):
      error = False
      if any([operation in item for operation in ['+', '-', ':']]):
        raise ValueError('Combining (mathematical) operators and symmetry' +
                         'labels is not supported')
      elif '*' in item:
        return [i for i,s in enumerate(self.get_sym()) if item.split('*.')[-1] in s]
      else:
        return [i for i,s in enumerate(self.get_sym()) if s == item]

    def parse_spin(item, all_alpha_beta):
      spindic = {0: 'all_alpha', 1: 'all_beta'}
      if isinstance(item, str):
        for i_s in range(2):
          if spindic[i_s] in item:
            all_alpha_beta[i_s] = True
            item = item.replace(spindic[i_s], '')
        if numpy.any(all_alpha_beta):
          return all_alpha_beta, item, None
        else:
          if 'alpha' in item:
            return all_alpha_beta, item.replace('alpha', ''), self.get_spin_index('alpha')
          elif 'beta' in item:
            return all_alpha_beta, item.replace('beta', ''), self.get_spin_index('beta')
          else:
            return all_alpha_beta, item, None
      elif isinstance(item, list):
        for i_s in range(2):
          if spindic[i_s] in item:
            all_alpha_beta[i_s] = True
            item = remove_from_list(item, spindic[i_s])
        if numpy.any(all_alpha_beta):
          return all_alpha_beta, item, None
        else:
          if 'alpha' in item:
            return all_alpha_beta, remove_from_list(item, 'alpha'), self.get_spin_index('alpha')
          elif 'beta' in item:
            return all_alpha_beta, remove_from_list(item, 'beta'), self.get_spin_index('beta')
          else:
            return all_alpha_beta, item, None

    
    regsplit = re.compile(r"[\s,;]")
    
    # We set these variables here for later reference
    all_alpha_beta = [False, False]

    if isinstance(fid_mo_list,str) and 'all_mo' in fid_mo_list.lower():
      all_alpha_beta, fid_mo_list, srec = parse_spin(fid_mo_list, all_alpha_beta)
      spinrestructions = [srec]
      mo_in_file_new = [[i for i in range(len(self.data))]]
      
    else:
      if isinstance(fid_mo_list,str) and not path.exists(fid_mo_list):
        if fid_mo_list == 'occupied':
          display('All occupied orbitals have been selected:')
          fid_mo_list = [map(str,(self.get_occ() > 1e-5).nonzero()[0])]
        elif fid_mo_list == 'unoccupied' or fid_mo_list == 'virtual':
          display('All unoccupied/virtual orbitals have been selected:')
          fid_mo_list = [map(str,(self.get_occ() <= 1e-5).nonzero()[0])]
        elif ',' in fid_mo_list:
          fid_mo_list = fid_mo_list.split(',')
        else:
          fid_mo_list = [fid_mo_list]
      if isinstance(fid_mo_list, list):
        for i in fid_mo_list:
          if not isinstance(i, list):
            if isinstance(i, int) or isinstance(i, float):
              i = str(int(i))
            i = re.sub(' +',' ', i)
            i = regsplit.split(i) if isinstance(i,str) else [i]
          mo_in_file.append(list(map(str,i)))
      else:
        try:
          fid=open(fid_mo_list,'r')
          flines = fid.readlines()
          fid.close()
          for line in flines:
            integer = line.replace(',',' ').split()
            mo_in_file.append(integer)
        except:
          raise IOError('The selected mo_list (%(m)s) is not valid!' %
                        {'m': fid_mo_list} + '\ne.g.\n\t1\t3\n\t2\t7\t9\n')
        fid_mo_list = mo_in_file

      # Print some information
      for i,j in enumerate(mo_in_file):
        display('\tLine %d: %s' % (i+1,', '.join(j)))
      display('')
      
      # Check if the molecular orbitals are specified by symmetry
      # (e.g. 1.1 in MOLPRO nomenclature) or
      # by the number in the input file (e.g. 1)
      mo_in_file_new = []
      spinrestructions = []
      for sublist in mo_in_file:
        all_alpha_beta, sublist, srec = parse_spin(sublist, all_alpha_beta)
        spinrestructions.append(srec)
        tmp = []
        for item in sublist:
          if '.' not in item and '*' not in item:
            if isinstance(item, str):
              tmp.extend(parse_nosym(item.lower()))
            else:
              tmp.extend(item)
          else:
            tmp.extend(parse_sym(item))
        mo_in_file_new.append(tmp)
    
    mo_in_file = []
    for isub, sublist in enumerate(mo_in_file_new):
      if numpy.any(all_alpha_beta):
        spindic = {0: 'alpha', 1: 'beta'}
        for i_s, all_spin in enumerate(all_alpha_beta):
          if all_spin:
            selected = []
            for isel in range(len(sublist)):
              if sublist[isel] in self.get_spin_index(spindic[i_s]):
                selected.append(isel)
      else:
        if spinrestructions[isub] is not None:
          selected = []
          for isel in range(len(sublist)):
            if sublist[isel] in spinrestructions[isub]:
              selected.append(isel)
        else:
          selected = range(len(sublist))
      mo_in_file.append([sublist[i] for i in selected])

    if flatten_input:
      mo_in_file = [[item for sublist in mo_in_file for item in sublist]]
      if sort_indices:
        mo_in_file = [numpy.sort([item for sublist in mo_in_file for item in sublist])]
    
    if len(mo_in_file) == 1 and flatten_input:
      mo_spec = MOClass([])
      mo_spec.spinpolarized = self.spinpolarized
      for index in mo_in_file[0]:
        mo_spec.append(self.data[index])
      mo_spec.selected_mo = mo_in_file[0]
      mo_spec.selection_string = str(fid_mo_list[0])
      # Update alpha_index and beta_index in newly created mo_spec
      mo_spec.alpha_index = require([i for i,j in enumerate(mo_in_file[0]) if j in self.alpha_index], dtype='i')
      mo_spec.beta_index = require([i for i,j in enumerate(mo_in_file[0]) if j in self.beta_index], dtype='i')

    else:
      mo_spec = []
      for i, sublist in enumerate(mo_in_file):
        mo_sub = MOClass([])
        mo_sub.spinpolarized = self.spinpolarized
        for index in sublist:
          mo_sub.append(self.data[index])
        mo_sub.selected_mo = sublist
        mo_sub.selection_string = str(fid_mo_list[i])
        # Update alpha_index and beta_index in newly created mo_spec
        mo_sub.alpha_index = require([i for i,j in enumerate(sublist) if j in self.alpha_index], dtype='i')
        mo_sub.beta_index = require([i for i,j in enumerate(sublist) if j in self.beta_index], dtype='i')
        mo_spec.append(mo_sub)

    # Print some information
    display('\nThe following orbitals will be considered...')
    for i, sublist in enumerate(mo_in_file):
      display('\tLine %d: %s' % (i+1, str(sublist)))

    display('')
    return mo_spec def molpro_mo_order_ci(occ_info,mo_spec,irreps=None,nIRREP=None,order_sym=False):
  '''Orders the molecular orbitals according to the occupation information 
  of a CI class. (For processing MCSCF MOLPRO output)  
  
  **Parameters:**
  
    occ_info : ci.info["occ_info"]
      Contains a dict with the number of orbitals in the different 
      IRREPS for "core", "closed", "active", and "external" orbitals.
    mo_spec : list of dictionaries 
      See :ref:`Central Variables` for details.
  
  **Returns:**
  
    closed,active,external : list of mo_spec
      Contains mo_spec assigned to the "closed", "active", and "external" 
      orbitals.
  '''
  if nIRREP is None: nIRREP = len(occ_info['closed']) 
  if irreps is None: 
    irreps = dict(zip(map(str,range(1,nIRREP+1)),range(nIRREP)))
  else:
    nIRREP = len(irreps)
    irreps = dict(zip(irreps,range(nIRREP)))
  
  mo_index = [[] for i in range(nIRREP)]
  for i,i_mo in enumerate(mo_spec):
    info = i_mo['sym'].split('.')
    mo_index[irreps[info[1]]].append((i,int(info[0])))
  
  mo_sym = [[] for i in range(nIRREP)]
  for i in range(nIRREP):
    if mo_index[i] == []:
      continue
    mo_index[i] = numpy.array(mo_index[i])
    j = numpy.argsort(mo_index[i],axis=0)[:,1]
    mo_index[i] = mo_index[i][j]
    mo_sym[i] = MOClass([mo_spec[j] for j,k in mo_index[i]])
    #mo_sym[i].update()
  if order_sym:
    return mo_sym
  
  closed = MOClass()
  active = MOClass()
  external = MOClass()
  for i in range(nIRREP):
    a = occ_info['core'][i] + occ_info['closed'][i]
    closed.extend(mo_sym[i][:a])
    b = a+occ_info['active'][i]
    active.extend(mo_sym[i][a:b])
    external.extend(mo_sym[i][b:])
  #closed.update()
  #active.update()
  #external.update()
  return closed,active,external def psi4_detci(fname,select_run=None,threshold=0.0,**kwargs):
  '''Reads PSI4 DETCI output. 
  
  **Parameters:**
  
    fname: str, file descriptor
      Specifies the filename for the input file.
      fname can also be used with a file descriptor instad of a filename.
    select_run : (list of) int
      Specifies the DETCI calculation (D E T C I) to be read. 
      For the selected DETCI calculation, all electronic states will be read.
    threshold : float, optional
      Specifies a read threshold for the CI coefficients.
  
  **Returns:**
  
    ci : list of CIinfo class instances 
      See :ref:`Central Variables` for details.
  
  ..#ATTENTION: Changed return value to list of CI classes
  '''

  display('\nReading data of DETCI calculation from PSI4...')
  count = 0
  numci = []
  # Go through the file line by line 

  if isinstance(fname, str):
    filename = fname
    fname = descriptor_from_file(filename, index=0, ci_descriptor=True)
  else:
    filename = fname.name

  from io import TextIOWrapper
  if isinstance(fname, TextIOWrapper):
    flines = fname.readlines()      # Read the WHOLE file into RAM
  else:
    magic = 'This is an Orbkit magic string'
    text = fname.read().decode("iso-8859-1").replace('\n','\n{}'.format(magic))
    flines = text.split(magic)
    flines.pop()

  for line in flines:
    if 'D E T C I' in line:
      count += 1
      numci.append(0)
    if '* ROOT' in line:
      numci[-1] += 1

  if count == 0:
    display('The input file %s does not contain any DETCI calculations!\n' % (filename) + 
            'It does not contain the keyword:\n\tD E T C I')
    raise IOError('Not a valid input file')
  else:
    string = ', '.join(map(str,numci)).rsplit(', ',1)
    string = ' and '.join(string) if len(string[0].split(',')) < 2 else ', and '.join(string)
    display('The input file %s contains' % (filename))
    display('%d DETCI calculation(s) with %s root(s)%s.'%(count, string,
                               ', respectively' if len(numci)>1 else ''))
    
    if select_run is None:
      select_run = numpy.arange(count)   
    if isinstance(select_run,int) and 0 <= select_run < count:
      select_run = [select_run]
      ci = [[]]
    elif isinstance(select_run,(list,numpy.ndarray)):
      ci = []
      for i in select_run:
        ci.append([])
        if not isinstance(i,(int, numpy.int, numpy.int64, numpy.intc)):
          raise IOError(str(i) + ' is not a valid selection for select_run')
    else: 
      raise IOError(str(select_run) + ' is a not valid selection for select_run')    
  select_run = numpy.array(select_run)
  
  display('\n\tYour selection (counting from zero): %s' % 
          ', '.join(map(str,select_run)))
  
  general_information = {'fileinfo': filename,
                         'read_threshold': threshold}
  
  ci_skip = 0
  count = 0
  count_runs = 0
  min_c = 1
  orbs_in_ci = 0
  start_reading = False
  
  occ_types = ['core','closed','active','external']
  synonyms = {'dropped docc': 'core',
              'frozen docc': 'core',
              'restricted docc': 'closed',
              'ras ': 'active',
              'active': 'active',
              'restricted uocc': 'external',
              'frozen uocc': 'external',
              'dropped uocc': 'external'  }
  irreps = []
  fist_active_mo = []
  index_active_mo = []
  occ_symbols = {'A': numpy.array([1,0]),
                 'B': numpy.array([0,1]),
                 'X': numpy.array([1,1])}

  flines = iter(flines)
  for line in flines:
    thisline = line.split()             # The current line split into segments
    
    if 'Running in ' in line and 'symmetry.' in line: 
      nIRREP = point_groups()[thisline[2].lower()]
      rhf_occ = numpy.zeros(nIRREP,dtype=numpy.intc)
    #--- RHF occupation
    elif 'Final Occupation by Irrep:' in line:
      if sys.version_info.major == 2:
        irreps = next(flines).split()
        line = next(flines)
      else:
        irreps = next(flines).split()
        line = next(flines)
      c_occ = line.replace(',','').split()[2:-1]
      for ii in range(len(c_occ)):
        rhf_occ[ii] = int(c_occ[ii])
    #--- A DETCI Calculation starts ---
    elif 'D E T C I' in line:
      occ_info = {}
      for i in occ_types:
        occ_info[i] = numpy.zeros(nIRREP,dtype=numpy.intc)
      info = {}
      num_roots = 0
      method = 'detci'
      i = numpy.argwhere(select_run == count_runs)
      if len(i) > 0:
        index_run = int(i)
        start_reading = True
      count_runs += 1
    elif start_reading:
      if 'NUM ROOTS     =' in line:
        num_roots = int(thisline[3])
      elif 'FCI          =' in line and line.endswith('yes'):
        method = 'fci'
      elif 'REF SYM       =' in line:
        info['sym'] = thisline[-1]
        if info['sym'] != 'auto':
          info['sym'] = irreps[int(info['sym'])]
      elif ' S     ' in line and ' =' in  line:
        info['spin'] = multiplicity()[int(2*float(thisline[2])+1)]
      elif 'Electrons    =' in line:
        info['nel'] = int(thisline[-1])
      elif 'ORBS IN CI   =' in line: 
        orbs_in_ci = int(thisline[-1])
      elif  any([i in line.lower() for i in synonyms.keys()]):
        for i in synonyms.keys():
          if i in line.lower():
            occ_info[synonyms[i]] += numpy.array(thisline[-nIRREP:],
                                                 dtype=numpy.intc)
      elif ('* ROOT' in line and 'CI total energy' in line) or 'CI Root' in line and 'energy' in line:
        if '*' in line:
          state = '%s.%s'%(thisline[2],info['sym'])
          energy = float(thisline[7])
        else:
          state = '%s.%s'%(int(thisline[2])+1,info['sym'])
          energy = float(thisline[-1])
        ci[index_run].append(CIinfo(method=method))
        ci[index_run][-1].info = copy(general_information) 
        ci[index_run][-1].info['fileinfo'] += '@%d' % index_run
        ci[index_run][-1].info['irreps'] = irreps
        ci[index_run][-1].info['state'] = state
        ci[index_run][-1].info['energy'] = energy
        ci[index_run][-1].info['spin'] = info['spin']
        ci[index_run][-1].info['nel'] = sum(rhf_occ)
        ci[index_run][-1].info['occ_info'] = occ_info
        closed = []
        active = {}
        c = 0
        for i in range(nIRREP):
          a = occ_info['core'][i] + occ_info['closed'][i]
          for b in range(occ_info['active'][i]):
            active['%d%s' % (a+b+1,irreps[i])] = c
            c += 1
      elif 'most important determinants' in line:
        num_det = int(thisline[1])
        ci[index_run][-1].coeffs = [] # numpy.zeros(num_det)
        ci[index_run][-1].occ = [] # numpy.zeros((numpy.sum(occ_info['active']),2),
                                        #dtype=numpy.intc)
        if sys.version_info.major == 2:
          next(flines)
        else:
          next(flines)
        def rm(s,w='*(,)'):
          for i in w:
            s = s.replace(i,' ')
          return s
        for i in range(num_det):  
          if sys.version_info.major == 2:
            thisline = rm(next(flines)).split()
          else:
            thisline = rm(next(flines)).split()
          if thisline == []:
            break
          min_c = min(min_c,abs(float(thisline[1])))
          if abs(float(thisline[1])) > threshold:
            ci[index_run][-1].coeffs.append(float(thisline[1]))
            ci[index_run][-1].occ.append(numpy.zeros((numpy.sum(occ_info['active']),2),
                                          dtype=numpy.intc))
            for j in thisline[4:]:
              ci[index_run][-1].occ[-1][active[j[:-1]]] = occ_symbols[j[-1]]
        ci[index_run][-1].coeffs = numpy.array(ci[index_run][-1].coeffs)
        ci[index_run][-1].occ = numpy.array(ci[index_run][-1].occ,dtype=numpy.intc)
      elif 'A good bug is a dead bug' in line:
        start_reading = False
  #--- Calculating norm of CI states
  display('\nIn total, %d states have been read.' % sum([len(i) for i in ci])) 
  display('Norm of the states:')
  for i in range(len(ci)):
    for j in range(len(ci[i])):
      norm = sum(ci[i][j].coeffs**2)
      display('\tState %s (%s):\tNorm = %0.8f (%d Coefficients)' % 
              (ci[i][j].info['state'],ci[i][j].info['spin'],
               norm,len(ci[i][j].coeffs)))
  display('')
  if threshold and min_c > threshold:
    display('\nInfo:'+
       '\n\tSmallest coefficient (|c|=%f) larger than the read threshold (%f).' 
       %(min_c,threshold) + 
       '\n\tUse `set num_dets_print -1` in the PSI4 input file to print ' +
       'all CI coefficients.\n')
  
  ci_new = []
  for i in ci:
    ci_new.extend(i)
  
  return ci_new def tmol_escf(filename,ene_unit='nm',read_props=False):
  '''Reads Turbomole TD-DFT output (``escf`` file). 

  **Parameters:**

  filename : str
    Specifies the filename for the input file.
  ene_unit : str
    Specifies in which units the excitation energies are read out.
    [nm,E_h,eV,cm-1]
  read_props : bool
    If True, reads the transition dipole moments and dominant molecular 
    orbital contributions.

  **Returns:**
  
  :if read_props: 
    - ex_ene,osc_str,tdm,dom_contrib
  :else:
    - ex_ene,osc_str
  
  ex_ene : numpy.ndarray, shape=(nstates)
    Contains the excitation energies for excited states.
  osc_str : numpy.ndarray, shape=(nstates)
    Contains the excitation energies for excited states.
  tdm : numpy.ndarray, shape=((nstates,) + 3)
    Contains the transition dipole moments from ground state to excited 
    states.
  dom_contrib : list, shape=((NMO,) + N)
    Contains the dominant molecular orbital contribution of excited states.
  '''

  fid    = open(filename,'r')      # Open the file
  flines = fid.readlines()         # Read the WHOLE file into RAM
  fid.close()                      # Close the file
  
  # Is this really a escf file? 
  if not '                                e s c f\n' in flines:
    raise IOError('The input file %s is no valid escf file!\n\nIt does'  % filename+
          ' not contain the keyword: e s c f\n')
  
  # Initialize variables
  props = {}
  props['E_0i'] = []
  props['dom_contrib'] = []
  props['f_0i'] = []
  props['mu_0i'] = []
  props['spin_i'] = []
  props['sym_i'] = []
  ncount = 1
  sec_flag = None
  
  # Select the excitation energy
  if ene_unit == 'nm':
    ene_str = 'Excitation energy / nm:'
  elif ene_unit == 'E_h':
    ene_str = 'Excitation energy:'
  elif ene_unit == 'eV':
    ene_str = 'Excitation energy / eV:'
  elif ene_unit == 'cm-1':
    ene_str = 'Excitation energy / cm^(-1):'
  else: 
    display('The units selection for the excitation energies is not valid.' +
            'The excitation energies are read out in units of Hartree.\n')
    ene_str = 'Excitation energy:'

  for il in range(len(flines)):
    line = flines[il]            # The current line as string
    thisline = line.split()      # The current line split into segments
    # Check the file for keywords 
    if ene_str in line:
      props['E_0i'].append(float(thisline[-1]))
    elif ('%s ' % ncount) in line and 'excitation' in line:
        props['sym_i'].append(thisline[2])
        props['spin_i'].append(thisline[1])
        ncount += 1
    elif 'Oscillator strength:' in line:
      sec_flag = 'osc'
      count = 1
    elif 'Dominant contributions:' in line:
      sec_flag = 'dcon'
      count = 2
      props['dom_contrib'].append([])
    elif 'Electric transition dipole moment (length rep.):' in line:
      sec_flag = 'mu_0i'
      count = 1
    else:
      if sec_flag == 'osc':
        if count == 0:
          props['f_0i'].append(float(thisline[-1]))
          sec_flag = None
        else:
          count -= 1
      elif sec_flag == 'mu_0i':
        if count == 0:
          tmp = numpy.zeros(3)
          for xyz in range(3):
            tmp[xyz] = float(flines[il+xyz].split()[1])
          props['mu_0i'].append(tmp)
          sec_flag = None
        else:
          count -= 1
      elif sec_flag == 'dcon':
        if count == 0:
          if len(thisline) == 7:
            tmp = []
            tmp.append(int(thisline[0]))
            tmp.append(int(thisline[3]))
            tmp.append(float(thisline[-1]))
            props['dom_contrib'][-1].append(tmp)
          else:
            sec_flag = None
        else:
          count -= 1
  for i in props.keys():
      props[i] = numpy.array(props[i])

  return props def tmol_tddft(fname,nmoocc=None,nforbs=0,select_state=None,threshold=0.0,
               bortho=False,**kwargs):
  '''Reads Turbomole TD-DFT output (``sing_a`` file). 
  
  Hint: Only implemented for singlet TD-DFT computations in C_1 symmetry, the
  output data file of which is usually called ``sing_a``.
  
  **Parameters:**
  
    fname: str, file descriptor
      Specifies the filename for the input file.
      fname can also be used with a file descriptor instad of a filename.
    nmoocc : int
      Specifies the number of non-frozen occupied orbitals.
    nforbs : int,optional
      Specifies the number of frozen orbitals.
    select_state : None or list of int, optional
      If not None, specifies the states to be read (0 corresponds to the ground 
      state), else read all electronic states. 
    threshold : float, optional
      Specifies a read threshold for the TD-DFT coefficients.
    bortho : bool
      If True, an orthonormalization of  the TD-DFT coefficients is perfomed
      with Gram-Schmidt.
  
  **Returns:**
  
    ci : list of CIinfo class instances
      See :ref:`Central Variables` for details.
  '''
  display('\nReading data of TD-DFT calculation from Turbomole...')
  
  nel = (nmoocc + nforbs) * 2
  lumo = (nmoocc + nforbs)
  states = []
  if isinstance(select_state,int): select_state = [select_state]

  if isinstance(fname, str):
    filename = fname
    fname = descriptor_from_file(filename, index=0, ci_descriptor=True)
  else:
    filename = fname.name

  start = False
  for i,line in enumerate(fname):
    if '$tensor space dimension' in line:
      tspace = int(line.split()[-1])
    if 'eigenvalue' in line:
      start = True
      l = line.split()
      states.append({'eigenvalue': float(l[-1].replace('D','E')),
                     'coeffs': []})
    elif start:
      for i in range((len(line)-1)//20):
        states[-1]['coeffs'].append(float(line[i*20:(i+1)*20].replace('D','E')))
  
  for i in range(len(states)):
    if len(states[i]['coeffs']) == tspace*2:
      states[i]['coeffs'] = numpy.array(states[i]['coeffs']).reshape((2,nmoocc,-1))
      states[i]['xia'] = states[i]['coeffs'][0]
    elif len(states[i]['coeffs']) == tspace:
      states[i]['xia'] = numpy.array(states[i]['coeffs']).reshape((nmoocc,-1))
    else:
      raise IOError('Shape of coefficient matrix does not match with tensor space.')
    del states[i]['coeffs']

  coeffs = numpy.zeros((len(states)+1,numpy.prod(states[1]['xia'].shape)+1))
  for i in range(len(states)):
    coeffs[i+1][1:] = states[i]['xia'].reshape((-1,))
  coeffs[0,0] = 1.0
  
  # Set up configurations
  ci = []
  for i in range(len(states)+1):
    if select_state is None or i in select_state:
      ci.append(CIinfo(method='tddft'))
      ci[-1].coeffs = coeffs[i]
      ci[-1].occ    = []
      ci[-1].info = {'state': str(i),
                     'energy': 0.0 if i==0 else numpy.abs(states[i-1]['eigenvalue']),
                     'fileinfo': filename,
                     'read_threshold': threshold,
                     'spin': 'Unknown',
                     'nel': nel}
      if not i:
        ci[-1].occ.append([-1,-1])
      else:
        ci[-1].occ.append([0,0])
        for jj in range(states[i-1]['xia'].shape[0]):
          for kk in range(states[i-1]['xia'].shape[1]):
            ci[-1].occ.append([jj+nforbs,kk+lumo])
      ci[-1].occ = numpy.array(ci[-1].occ,dtype=numpy.intc)

  # Gram-Schmidt
  display('Orthonormalizing the TD-DFT coefficients with Gram-Schmidt...\n')
  ci = orthonorm(ci,reorth=False)

  if bortho:
    for c in ci:
      c.apply_threshold(threshold,keep_length=True)
    ci = orthonorm(ci,reorth=True)
  else:
    for c in ci:
      c.apply_threshold(threshold,keep_length=False)
  
  for c in ci:
      c.apply_threshold(0.0,keep_length=False)
    
  #--- Calculating norm of CI states
  display('\nIn total, %d states have been read.' % len(ci)) 
  display('Norm of the states:')
  for i in ci:
    display(str(i))
  
  return ci def gamess_cis(filename,select_state=None,threshold=0.0,**kwargs):
  '''Reads GAMESS-US CIS output. 
  
  **Parameters:**
  
    filename : str
      Specifies the filename for the input file.
    select_state : None or list of int, optional
      If not None, specifies the states to be read (0 corresponds to the ground 
      state), else read all electronic states.
    threshold : float, optional
      Specifies a read threshold for the CI coefficients.
  
  **Returns:**
  
    ci : list of CIinfo class instances
      See :ref:`Central Variables` for details.
  '''
  display('\nReading data of CIS calculation from GAMESS-US...')
  # Initialize variables
  ci = []
  ci_flag = False
  prttol = False
  init_state = False
  rhfspin = 0
  min_c = -1
  
  if isinstance(select_state,int): select_state = [select_state]
  with open(filename) as fileobject:
    for line in fileobject:
      thisline = line.split()             # The current line split into segments
      #--- Check the file for keywords ---
      # Initialize Hartree-Fock ground state
      if 'NUMBER OF ELECTRONS' in line and "=" in line:
        nel = int(thisline[-1])
      elif 'SPIN MULTIPLICITY' in line:
        rhfspin = int(thisline[-1])
      elif ' FINAL RHF ENERGY IS' in line and (select_state is None or 0 in select_state):
          ci.append(CIinfo(method='cis'))
          ci[-1].info   = []
          ci[-1].coeffs = []
          ci[-1].occ    = []
          ci[-1].occ.append([0,0])
          ci[-1].coeffs.append(1.0)
          ci[-1].info = {'state': '0',
                         'energy': float(thisline[4]),
                         'fileinfo': filename,
                         'read_threshold': threshold,
                         'spin': multiplicity()[rhfspin],
                         'nel': nel}
      # Printing parameter
      elif ' PRINTING CIS COEFFICIENTS LARGER THAN' in line:
        min_c = float(thisline[-1])
      # Initialize new excited state
      elif ' EXCITED STATE ' in line and 'ENERGY=' and 'SPACE SYM' in line:
        if select_state is None or int(thisline[2]) in select_state:
          init_state = True
          cis_skip = 6
          ci.append(CIinfo(method='cis'))
          ci[-1].info   = []
          ci[-1].coeffs = []
          ci[-1].occ    = []
          ci[-1].info = {'state': thisline[2],
                         'energy': float(thisline[4]),
                         'fileinfo': filename,
                         'read_threshold': threshold,
                         'spin': multiplicity()[int(2*float(thisline[7])+1)],
                         'nel': nel}
      if init_state == True:
        if not cis_skip:
          if '----------------------------------------------' in line:
            init_state = False
          else:
            if abs(float(thisline[2])) > threshold:
              ci[-1].occ.append(thisline[:2])
              ci[-1].coeffs.append(thisline[2])
        elif cis_skip:
          cis_skip -= 1
          
  #--- Calculating norm of CI states
  display('\nIn total, %d states have been read.' % len(ci)) 
  display('Norm of the states:')
  for i in range(len(ci)):
    j = numpy.array(ci[i].coeffs,dtype=float)
    norm = numpy.sum(j**2)
    ci[i].coeffs = j
    # Write Norm to log-file
    display('\tState %s:\tNorm = %0.8f (%d Coefficients)' % (ci[i].info['state'],norm, len(ci[i].coeffs)))
    # Transform to numpy arrays
    ci[i].occ = numpy.array([s for s in ci[i].occ],dtype=numpy.intc)-1
  display('')
  if min_c > threshold:
    display('\nInfo:'+
       '\n\tSmallest coefficient (|c|=%f) is larger than the read threshold (%f).' 
       %(min_c,threshold) + 
       '\n\tUse `PRTTOL=0.0` in the `$CIS` input card to print ' +
       'all CI coefficients.\n')
  
  return ci def gamess_tddft(fname,select_state=None,threshold=0.0,**kwargs):
  '''Reads GAMESS-US TDDFT output. 
  
  **Parameters:**
  
    fname: str, file descriptor
      Specifies the filename for the input file.
      fname can also be used with a file descriptor instad of a filename.
    select_state : None or list of int, optional
      If not None, specifies the states to be read (0 corresponds to the ground 
      state), else read all electronic states.
    threshold : float, optional
      Specifies a read threshold for the CI coefficients.
  
  **Returns:**
  
    ci : list of CIinfo class instances
      See :ref:`Central Variables` for details.
  '''
  display('\nReading data of TDDFT calculation from GAMESS-US...')
  # Initialize variables
  ci = []
  ci_flag = False
  prttol = False
  init_state = False
  rhfspin = 0
  spin = 'Unknown'
  
  if isinstance(select_state,int): select_state = [select_state]

  if isinstance(fname, str):
    filename = fname
    fname = descriptor_from_file(filename, index=0, ci_descriptor=True)
  else:
    filename = fname.name

  for line in fname:
    thisline = line.split()             # The current line split into segments
    #--- Check the file for keywords ---
    # Initialize Hartree-Fock ground state
    if 'NUMBER OF ELECTRONS' in line:
      nel = int(thisline[-1])
    elif 'SPIN MULTIPLICITY' in line:
      rhfspin = int(thisline[-1])
    elif 'SINGLET EXCITATIONS' in line:
      spin = 'Singlet'
    #elif ' FINAL RHF ENERGY IS' in line and (select_state is None or 0 in select_state):
    elif ' FINAL' in line and ' ENERGY IS' in line and (select_state is None or 0 in select_state):
        ci.append(CIinfo(method='tddft'))
        ci[-1].info   = []
        ci[-1].coeffs = []
        ci[-1].occ    = []
        ci[-1].occ.append([0,0])
        ci[-1].coeffs.append(1.0)
        ci[-1].info = {'state': '0',
                       'energy': float(thisline[4]),
                       'fileinfo': filename,
                       'read_threshold': threshold,
                       'spin': spin,
                       'nel': nel}
    # Initialize new excited state
    elif 'STATE #' in line and 'ENERGY =' in line:
      if select_state is None or int(thisline[2]) in select_state:
        init_state = True
        tddft_skip = 8
        ci.append(CIinfo(method='cis'))
        ci[-1].info   = []
        ci[-1].coeffs = []
        ci[-1].occ    = []
        ci[-1].info = {'state': thisline[2],
                       'energy': float(thisline[-2])*ev_to_ha + ci[0].info['energy'],
                       'fileinfo': filename,
                       'read_threshold': threshold,
                       'spin': 'Unknown',
                       'nel': nel}
    if init_state == True and line != '\n' and 'WARNING:' not in line:
      if not tddft_skip:
        if 'NON-ABELIAN' in line or 'SUMMARY' in line or 'SYMMETRY' in line or 'STATE #' in line:
          init_state = False
        else:
          if abs(float(thisline[2])) > threshold:
            ci[-1].occ.append(thisline[:2])
            ci[-1].coeffs.append(thisline[2])
      elif tddft_skip:
        tddft_skip -= 1
          
  fname.close()
          
  #--- Calculating norm of CI states
  display('\nIn total, %d states have been read.' % len(ci)) 
  display('Norm of the states:')
  for i in range(len(ci)):
    j = numpy.array(ci[i].coeffs,dtype=float)
    norm = numpy.sum(j**2)
    ci[i].coeffs = j
    # Write Norm to log-file
    display('\tState %s:\tNorm = %0.8f (%d Coefficients)' % (ci[i].info['state'],norm, len(ci[i].coeffs)))
    # Transform to numpy arrays
    ci[i].occ = numpy.array([s for s in ci[i].occ],dtype=numpy.intc)-1
  
  return ci def molpro_mcscf(fname,select_run=0,threshold=0.0,**kwargs):
  '''Reads MOLPRO MCSCF output. 
  
  **Parameters:**
  
    fname: str, file descriptor
      Specifies the filename for the input file.
      fname can also be used with a file descriptor instad of a filename.
    select_run : (list of) int
      Specifies the MCSCF calculation (1PROGRAM * MULTI) to be read. 
      For the selected MCSCF calculation, all electronic states will be read.
    threshold : float, optional
      Specifies a read threshold for the CI coefficients.
  
  **Returns:**
  
    ci : list of CIinfo class instances 
      See :ref:`Central Variables` for details.
  

    ..ATTENTION: Changed return value to list of CI classes
  '''
  display('\nReading data of MCSCF calculation from MOLPRO...')
  method = 'mcscf'
  available = {'mcscf': 'MULTI'}
  single_run_selected = isinstance(select_run,int)
  count = 0
  numci = []

  if isinstance(fname, str):
    filename = fname
    fname = descriptor_from_file(filename, index=0, ci_descriptor=True)
  else:
    filename = fname.name

  from io import TextIOWrapper
  if isinstance(fname, TextIOWrapper):
    flines = fname.readlines()      # Read the WHOLE file into RAM
  else:
    magic = 'This is an Orbkit magic string'
    text = fname.read().decode("iso-8859-1").replace('\n','\n{}'.format(magic))
    flines = text.split(magic)
    flines.pop()

  # Go through the file line by line 
  for line in flines:
    if '1PROGRAM * %s' % available[method] in line:
      count += 1
      numci.append(0)
    if '!%s state' % method in line.lower() and 'energy' in line.lower():
      numci[-1] += 1
          
  if count == 0:
    display('The input file %s does not contain any DETCI calculations!\n' % (filename) + 
            'It does not contain the keyword:\n\t1PROGRAM * MULTI')
    raise IOError('Not a valid input file')
  else:
    string = ', '.join(map(str,numci)).rsplit(', ',1)
    string = ' and '.join(string) if len(string[0].split(',')) < 2 else ', and '.join(string)
    display('The input file %s contains' % (filename))
    display('%d MCSCF calculation(s) with %s root(s)%s.'%(count,string,
                               ', respectively' if len(numci)>1 else ''))
    
    if select_run is None:
      select_run = numpy.arange(count)   
    if isinstance(select_run,int) and 0 <= select_run < count:
      select_run = [select_run]
      ci = [[]]
    elif isinstance(select_run,(list,numpy.ndarray)):
      ci = []
      for i in select_run:
        ci.append([])
        if not isinstance(i,int):
          raise IOError(str(i) + ' is not a valid selection for select_run')
    else: 
      raise IOError(str(select_run) + ' is a not valid selection for select_run')    
  select_run = numpy.array(select_run)
  
  display('\tYour selection (counting from zero): %s' % 
          ', '.join(map(str,select_run)))
  
  general_information = {'fileinfo': filename,
                         'read_threshold': threshold}
  
  ci_skip = 0
  count = 0
  count_runs = 0
  min_c = 1
  start_reading = False
  sec_flag = False

  info_flag = False
  info_sep = False
  info_split = False
    
  occ_types = ['core','closed','active','external']
  
  flines = iter(flines)
  for line in flines:
    thisline = line.split()             # The current line split into segments
    
    if '1PROGRAM *' in line:
      start_reading = False
    #--- Number of IRREPs
    if 'Point group' in line or '_PGROUP' in line: 
      nIRREP = point_groups()[thisline[-1].lower()]
      rhf_occ = numpy.zeros(nIRREP,dtype=numpy.intc)
    #--- RHF occupation
    elif 'Final occupancy:' in line:
      c_occ = line.split()[2:]
      for ii in range(len(c_occ)):
        rhf_occ[ii] = int(c_occ[ii])
    #--- A MCSCF Calculation starts ---
    elif '1PROGRAM * %s' % available[method] in line:
      occ_info = {}
      for i in occ_types:
        occ_info[i] = numpy.zeros(nIRREP,dtype=numpy.intc)
      state_info = []
      i = numpy.argwhere(select_run == count_runs)
      if len(i) > 0:
        index_run = int(i)
        start_reading = True
        count = 0
        old = 0
      count_runs += 1
    elif start_reading:
      #--- Active space ---
      if 'Number of ' in line and 'orbitals:' in line:          
        line = line.replace('(','').replace(')','').replace('-shell','')   
        c_occ = numpy.array(line.split()[-nIRREP:],dtype=numpy.intc)
        occ_info[line.split()[2]] += c_occ
      elif 'State symmetry' in line:
        next(flines)
        thisline = next(flines)
        if 'State symmetry' in thisline:
          next(flines)
          thisline = next(flines)
        thisline = thisline.replace('=',' ').split()
        data = {'nel': thisline[3],
                'spin': thisline[6],
                'sym': thisline[9]}
        thisline = next(flines).split()
        state_info.extend([data for i in range(int(thisline[-1]))])
      elif '!%s state' % method in line.lower() and 'energy' in line.lower():
        ci[index_run].append(CIinfo(method=method))
        info = state_info[count]
        thisline = line.lower().replace('state','state ').split()
        ci[index_run][count].info = copy(general_information) 
        ci[index_run][count].info['fileinfo'] += '@%d' % index_run
        ci[index_run][count].info['state'] = thisline[2]
        ci[index_run][count].info['energy'] = float(thisline[4])
        ci[index_run][count].info['spin'] = info['spin']
        ci[index_run][count].info['nel'] = info['nel']
        ci[index_run][count].info['occ_info'] = occ_info
        count += 1
      elif 'CI vector' in line:
        sec_flag = 'mcscf'
        info_split = '     '
        ci_skip = 3
        info = thisline[-1]
        count = old
        first = True
      if not ci_skip:
        if line == '\n' or '/EOF' in line:
          sec_flag = False
        elif sec_flag != False:
          split_line = list(filter(None, line.split(info_split)))
          if len(split_line) > 1:
            occupation = numpy.zeros((numpy.sum(occ_info['active']),2),dtype=numpy.intc)
            for i,j in enumerate(split_line[0].replace(' ','')):
              if j == '2': 
                occupation[i,:] = 1
              elif j == 'a':
                occupation[i,0] = 1
              elif j == 'b':
                occupation[i,1] = 1
            c0 = 0
          coeffs = split_line[-1].split()
          for i,j in enumerate(coeffs):
            if first:
              old += 1
            min_c = min(min_c,abs(float(j)))
            if abs(float(j)) > threshold:
              ci[index_run][count + c0 + i].coeffs.append(float(j))
              ci[index_run][count + c0 + i].occ.append(occupation)
          c0 += len(coeffs)
          first = False
      else:
        ci_skip -= 1 
  
  #--- Calculating norm of CI states
  display('\nIn total, %d states have been read.' % sum([len(i) for i in ci])) 
  display('Norm of the states:')
  for i in range(len(ci)):
    for j in range(len(ci[i])):
      ci[i][j].coeffs = numpy.array(ci[i][j].coeffs)
      ci[i][j].occ = numpy.array(ci[i][j].occ,dtype=numpy.intc)
      norm = sum(ci[i][j].coeffs**2)
        # Write Norm to log-file
      display('\tState %s (%s):\tNorm = %0.8f (%d Coefficients)' % 
              (ci[i][j].info['state'],ci[i][j].info['spin'],
               norm,len(ci[i][j].coeffs)))
  display('')
  if min_c > threshold:
    display('\nInfo:'+
       '\n\tSmallest coefficient (|c|=%f) larger than the read threshold (%f).' 
       %(min_c,threshold) + 
       '\n\tUse `gthresh, printci=0.0` in the MOLPRO input file to print ' +
       'all CI coefficients.\n')
  
  #if single_run_selected and len(ci) == 1:
    #ci = ci[0]
  ci_new = []
  for i in ci:
    ci_new.extend(i)
  
  return ci_new def gaussian_tddft(fname,select_state=None,threshold=0.0,**kwargs):
  '''Reads Gaussian16 TDDFT output. 
  
  **Parameters:**
  
    fname: str, file descriptor
      Specifies the filename for the input file.
      fname can also be used with a file descriptor instad of a filename.
    select_state : None or list of int, optional
      If not None, specifies the states to be read (0 corresponds to the ground 
      state), else read all electronic states.
    threshold : float, optional
      Specifies a read threshold for the CI coefficients.
  
  **Returns:**
  
    ci : list of CIinfo class instances
      See :ref:`Central Variables` for details.
  '''
  display('\nReading data of TDDFT calculation from Gaussian...')
  # Initialize variables
  ci = []
  ci_flag = False
  prttol = False
  init_state = False
  rhfspin = 0
  spin = 'Unknown'
  nel = 0
  deex = []
  
  if isinstance(select_state,int): select_state = [select_state]

  if isinstance(fname, str):
    filename = fname
    fname = descriptor_from_file(filename, index=0, ci_descriptor=True)
  else:
    filename = fname.name

  for line in fname:
    thisline = line.split()             # The current line split into segments
    #--- Check the file for keywords ---
    # Initialize Hartree-Fock ground state
    if ' SCF Done:  E(' in line:
        ci.append(CIinfo(method='tddft'))
        ci[-1].info   = []
        ci[-1].coeffs = []
        ci[-1].occ    = []
        ci[-1].occ.append([0,0])
        ci[-1].coeffs.append(1.0)
        ci[-1].info = {'state': '0',
                       'energy': float(thisline[4]),
		       'energy_nm': 0.0,
                       'fileinfo': filename,
                       'read_threshold': threshold,
                       'spin': spin,
		       'f_0i': 0.0}
    # Initialize new excited state
    elif ' Excited State' in line and 'eV' in line and 'nm' in line:
      if select_state is None or int(thisline[2].replace(':',' ')) in select_state:
        init_state = True
        tddft_skip = 1
        ci.append(CIinfo(method='tddft'))
        ci[-1].info   = []
        ci[-1].coeffs = []
        ci[-1].occ    = []
        ci[-1].info = {'state': thisline[2][:-1],
                       'energy': float(thisline[-6])*ev_to_ha + ci[0].info['energy'],
		       'energy_nm': float(thisline[-4]),
                       'fileinfo': filename,
                       'read_threshold': threshold,
                       'spin': thisline[3].split('-')[0],
		       'f_0i': float(thisline[8].replace('=',' ').split()[-1])}
        deex.append([])
    if init_state == True:
      if not tddft_skip:
        if thisline == [] or '->' not in line and '<-' not in line:
          init_state = False
        else:
          if '->' in line:
            thisline = line.replace('->','-> ').split()
            if abs(float(thisline[-1])) > threshold:
              tmp_occ = [thisline[0],thisline[2]]
              ci[-1].occ.append(tmp_occ)
              ci[-1].coeffs.append(float(thisline[-1])*numpy.sqrt(2))
          elif '<-' in line:
            deex[-1].append(float(thisline[-1])*numpy.sqrt(2))
      elif tddft_skip:
        tddft_skip -= 1

  fname.close()
  deex = numpy.array(deex)

  #--- Calculating norm of CI states
  display('\nIn total, %d states have been read.' % len(ci)) 
  display('Norm of the states:')
  for i in range(len(ci)):
    j = numpy.array(ci[i].coeffs,dtype=float)
    norm = numpy.sum(j**2)
    ci[i].coeffs = j
    # Write Norm to log-file
    display('\tState %s:\tNorm = %0.8f (%d Coefficients)' % (ci[i].info['state'],norm, len(ci[i].coeffs)))
    # Transform to numpy arrays
    ci[i].occ = numpy.array([s for s in ci[i].occ],dtype=numpy.intc)-1
  
  return ci
import numpy
from copy import copy
from orbkit.display import display
from orbkit.qcinfo import QCinfo, CIinfo
from orbkit.orbitals import MOClass
from orbkit.units import ev_to_ha
from .tools import molpro_mo_order_ci
from .psi4 import psi4_detci
from .tmol import tmol_escf, tmol_tddft
from .gamess import gamess_cis, gamess_tddft
from .molpro import molpro_mcscf
from .gaussian import gaussian_tddft



def main_ci_read(qc,fname,itype='psi4_detci',threshold=0.0,
                 select=None,nforbs=0,bortho=False,
                 **kwargs):
  '''Reads determinant CI calculation. 
  
  Supported input files (``itype``):
  
  itype           QC-Program  Level of Theory
  ==============  ==========  =====================
  'psi4_detci'    PSI4        All CI calculations
  'gamess_cis'    GAMESS-US   CIS
  'gamess_tddft'  GAMESS-US   TD-DFT
  'gauss16_tddft' Gaussian16  TD-DFT
  'tmol_tddft'    TURBOMOLE   TD-DFT
  'molpro_mcscf'  MOLPRO      MCSCF
  ==============  ==========  =====================
  
  **Parameters:**
  
    qc : class QCinfo
      See :ref:`Central Variables` for details.
    fname: str, file descriptor
      Specifies the filename for the input file.
      fname can also be used with a file descriptor instad of a filename.
    itype : str, choices={'tar', 'psi4_detci', 'gamess_cis', 'tmol_tddft', 'molpro_mcscf'}
      Specifies the type of the input file.
    threshold : float, optional
      Specifies a read threshold for the CI coefficients.
    select : None or (list of) int (counting from zero), optional
      if None, reads all states, else ...
      | TURBOMOLE & GAMESS-US: Specifies the STATE to be read (0 -> ground state)
      | PSI4 & MOLPRO: Specifies the CALCULATION to be read
    nforbs : int, optional, TURBOMOLE-specific
      Specifies the number of frozen orbitals.
    bortho : bool, optional, TURBOMOLE-specific
      If True, an orthonormalization of  the TD-DFT coefficients is perfomed
      with Gram-Schmidt.
  
  **Returns:**
  
    qc : class QCinfo
      New instance of QCinfo. See :ref:`Central Variables` for details.
      qc.mo_spec is possibly reordered.
    ci : list of CIinfo class instances 
      See :ref:`Central Variables` for details.
  '''

  if isinstance(fname, str):
    filename = fname
  else:
    filename = fname.name


  display('Opened \n\t%s\n' % filename)  