new_contents: #!/usr/bin/env python
  # vim: set fileencoding=utf-8 ts=8 sw=4 tw=0 :
  
  """
  Convert the discrete CP2K PDOS points to a smoothed curved using convoluted gaussians
  """
  
  # Copyright (c) 2017 Tiziano Müller <tiziano.mueller@chem.uzh.ch>,
  # based on a Fortran tool written by Marcella Iannuzzi <marcella.iannuzzi@chem.uzh.ch>
  #
  # This program is free software: you can redistribute it and/or modify
  # it under the terms of the GNU General Public License as published by
  # the Free Software Foundation, either version 3 of the License, or
  # (at your option) any later version.
  #
  # This program is distributed in the hope that it will be useful,
  # but WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  # GNU General Public License for more details.
  #
  # You should have received a copy of the GNU General Public License
  # along with this program. If not, see <http://www.gnu.org/licenses/>.
  
  
  from __future__ import print_function
  
  import sys
  import re
  import argparse
  
  import numpy as np
  
  
  HEADER_MATCH = re.compile(
      r'\# Projected DOS for atomic kind C at iteration step i = \d+, E\(Fermi\) = [ \t]* (?P<Efermi>[^\t ]+) a\.u\.')
  
  
  if __name__ == '__main__':
      parser = argparse.ArgumentParser(description=__doc__)
      parser.add_argument('pdosfilename', metavar='pdos-file', type=str,
                          help="the PDOS file generated by CP2K")
      parser.add_argument('natoms', metavar='number-of-atoms', type=int,
                          help="the number of atoms in the structure")
      parser.add_argument('--sigma', '-s', type=float, default=0.02,
                          help="sigma for the gaussian distribution")
      parser.add_argument('--de', '-d', type=float, default=0.001,
                          help="integration step size")
      parser.add_argument('--total-sum', action='store_true',
                          help="calculate and print the total sum for each orbital")
      args = parser.parse_args()
  
      with open(args.pdosfilename, 'r') as fhandle:
          header = HEADER_MATCH.match(fhandle.readline().rstrip())
          Efermi = float(header.group('Efermi'))  # will raise an error if not found
          data = np.loadtxt(fhandle)  # load the rest directly with numpy
  
      npnts, ncols = data.shape
      ncols -= 3  # ignore energy, occupation and #MO for mesh
  
      emin = min(data[:, 1]) - 0.2 # energies are in the second column
      emax = max(data[:, 1]) + 0.2
      nmesh = int((emax-emin)/args.de)+1
  
      # printing to stderr makes it possible to simply redirect the stdout to a file
      print("# of lines:       {:14d}".format(npnts), file=sys.stderr)
      print("Integration step: {:14.3f}".format(args.de), file=sys.stderr)
      print("Sigma:            {:14.3f}".format(args.sigma), file=sys.stderr)
      print("Minimum energy:   {:14.3f}".format(emin), file=sys.stderr)
      print("Maximum energy:   {:14.3f}".format(emax), file=sys.stderr)
      print("# of mesh points: {:14d}".format(nmesh), file=sys.stderr)
  
      # reproduce exactly the previous Fortran-based code
      xmesh = np.array([emin+i*args.de for i in range(1, nmesh+1)])
      ymesh = np.zeros((nmesh, ncols))
  
      fact = args.de/(args.sigma*np.sqrt(2.0*np.pi))
      for mpnt in range(nmesh):
          func = np.exp(-(xmesh[mpnt]-data[:, 1])**2/(2.0*args.sigma**2))*fact
          ymesh[mpnt, :] = func.dot(data[:, 3:])
  
      if args.total_sum:
          finalsum = np.sum(ymesh, 0)*args.de
          print("Sum over all meshpoints, per orbital:", file=sys.stderr)
          print(("{:16.8f}"*ncols).format(*finalsum), file=sys.stderr)
  
      xmesh -= Efermi  # put the Fermi energy at 0
      xmesh *= 27.211384  # convert to eV
      ymesh /= args.natoms  # normalize
      for mpnt in range(nmesh):
          print(("{:16.8f}" + " {:16.8f}"*ncols).format(xmesh[mpnt], *ymesh[mpnt, :]))
  
  #  vim: set ts=4 sw=4 tw=0 :
  
subject: convert CP2K PDOS output to gaussian-smoothed plots