def grid_init(is_vector=False, force=False):
def get_grid(start='\t'):
def tolist():
def todict():
def get_shape():
def set_grid(xnew,ynew,znew,is_vector):
def set_boundaries(is_regular,Nx=None,Ny=None,Nz=None):
def get_bbox():
def grid2vector():
def vector2grid(Nx,Ny,Nz):
def matrix_grid2vector(matrix): 
def matrix_vector2grid(matrix,Nx=None,Ny=None,Nz=None): 
def mv2g(**kwargs):
def grid_sym_op(symop):
def grid_translate(dx,dy,dz):
def rot(ang,axis):
def reflect(plane):
def inversion():
def sph2cart_vector(r,theta,phi):
def cyl2cart_vector(r,phi,zed):
def random_grid(geo_spec,N=1e6,scale=0.5):
def read(filename, comment='#'):
  def check(i, is_vector):
def adjust_to_geo(qc,extend=5.0,step=0.1):
def check_atom_select(atom,geo_info,geo_spec,interactive=True,
                      display=sys.stdout.write):
def center_grid(ac,display=sys.stdout.write):
def reset_grid():
  N_ = [len(x),len(y),len(z)]
    N_ = [Nx,Ny,Nz]
  P=[numpy.zeros((3,1)), numpy.reshape(ac,(3,1))]
  N_   = [len(grid[0]), len(grid[1]), len(grid[2])]
  N_   = [ 101,  101,  101]
N_   = [ 101,  101,  101]   #: Specifies the number of grid points (regular grid). def get_options():
def init_parser():
  def default_if_called(option, opt, value, parser, default=1e4):
def raise_error(string,error=IOError):
def print_message(string):
def check_options(error=raise_error,display=print_message,
                  interactive=False,info=True,check_io=True):
  def check_mo(attr):
def check_if_exists(fid, what='',error=IOError,display=sys.stdout.write,
                interactive=False):
def check_grid_output_compatibilty(error=raise_error):  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() def amira_creator(data,filename):
  '''Creates a ZIBAmira mesh file. (plain text)
  
  **Parameters:**
  
  data : numpy.ndarray, shape=N or (3,)+N
    Contains the output data.
  filename : str
    Contains the base name of the output file.
  '''
  if data.ndim == 3:
    typename = 'double'
    formatter = '{0:.15e}\n'
    data = numpy.swapaxes(data,0,2).reshape((-1,1)) #ravel(order='F')[:,numpy.newaxis]
  elif data.ndim == 4:
    typename = 'double[3]'
    data = numpy.array([j.ravel(order='F') for j in data])
    formatter = '{0:.15e} {1:.15e} {2:.15e}\n'
  else:
    raise IOError("amira_creator only supports 3D or 4D data.")
  
  N = tuple(grid.N_)
  
  # Open an empty file 
  filename += '.am' if not filename.endswith('.am') else ''
  fid = open(filename,'w')
  # Write Header 
  fid.write('# AmiraMesh 3D ASCII 2.0\n\n\n')
  fid.write('define Lattice %d %d %d\n' % N)
  fid.write('Parameters {\n')
  fid.write('    Content "%dx%dx%d %s, uniform coordinates",\n' % 
            (N + (typename,)))
  fid.write('    BoundingBox %(xmin)f %(xmax)f %(ymin)f %(ymax)f %(zmin)f %(zmax)f,\n' %
            {'xmin': grid.min_[0],'xmax': grid.max_[0],
             'ymin': grid.min_[1],'ymax': grid.max_[1],
             'zmin': grid.min_[2],'zmax': grid.max_[2]})
  fid.write('    CoordType "uniform"\n}\n\n')
  fid.write('Lattice { %s Data } @1\n'%typename)
  fid.write('# Data section follows\n@1\n')
  for i in data:
    fid.write(formatter.format(*i))
  #for tt in range(len(grid.z)):
    #for ss in range(len(grid.y)):
      #for rr in range(len(grid.x)): 
        #string += '%g\n' % rho[rr,ss,tt]
  fid.close() def hx_network_creator(rho,filename):
  '''Creates a ZIBAmira hx-network file including a colormap file (.cmap)
  adjusted to the density for the easy depiction of the density.
  '''
  from .hx_network_draft import hx_network
  filename += '.hx' if not filename.endswith('.hx') else ''
    
  # Create a .cmap colormap file using the default values 
  display('\tCreating ZIBAmira colormap file...\n\t\t' + filename.replace('.hx','.cmap'))
  
  assert (rho.shape != tuple(grid.N_)), 'The grid does not fit the data.'
  
  colormap_creator(rho,filename.replace('.hx','.cmap'))
  
  # Create a .hx network file based on the file orbkit.hx_network_draft.py 
  # Open an empty file
  
  fid = open(filename,'w')
  
  # Copy the content of the draft file and replace the keywords 
  fid.write(hx_network.replace("FILENAME",os.path.splitext(os.path.basename(filename))[0])) 
  
  # Close the file 
  fid.close()   def cube_creator(data,filename,geo_info,geo_spec,comments='',labels=None,**kwargs):
  '''Creates a plain text Gaussian cube file. 
  
  **Parameters:**
  
  data : numpy.ndarray, shape=N
    Contains the output data.
  filename : str
    Contains the base name of the output file.
  geo_info, geo_spec : 
    See :ref:`Central Variables` for details.
  comments : str, optional
    Specifies the second (comment) line of the cube file.  
  '''
  # Shape shall be (Ndrv,Ndata,Nx,Ny,Nz) or (Ndrv,Ndata,Nxyz)
  data = numpy.array(data)
  dims = 3
  ndim = data.ndim
  shape = data.shape
    
  if data.ndim < dims:
    raise AssertionError('data.ndim < ndim of grid')
  elif data.ndim == dims: # 3d data set
    data = data[numpy.newaxis]
  elif data.ndim > dims + 1:
    raise AssertionError('data.ndim > (ndim of grid) +2')
  
  if labels is not None:
    if labels == True or labels == 'auto':
      labels = list(range(len(data)))
    assert len(labels) == len(data)
    try: 
      labels = map(int,labels)
    except ValueError:
      raise AssertionError('labels has to be list of integers.')

  assert data.shape[1:] == tuple(grid.N_), 'The grid does not fit the data.'
  
  if not any([filename.endswith(ext)
              for ext in ['cube', 'cb', 'cube.gz', 'cb.gz']]):
    filename += '.cube'
  
  # Write the type and the position of the atoms in the header
  string = 'orbkit calculation\n'
  string += ' %(f)s\n'  % {'f': comments}
  # How many atoms 
  string += ('%(at)d' % {'at': (-1)**(labels is not None)*len(geo_info)}).rjust(5)
  # Minima
  for ii in range(3):
    string += ('%(min)0.6f' % {'min': grid.min_[ii]}).rjust(12)
  # Number of data points per grid point
  if len(data) > 1:
    string += ('%d' % len(data)).rjust(12)
  for ii in range(3):
    string += '\n'
    string += ('%(N)d' % {'N': grid.N_[ii]}).rjust(5)
    for jj in range(3):
      if jj == ii: 
        string += ('%(dr)0.6f' % {'dr': grid.delta_[ii]}).rjust(12)
      else:
        string += ('%(dr)0.6f' % {'dr': 0}).rjust(12)
  
  string += '\n'
  for ii in range(len(geo_info)):
    string += ('%(N)d' % {'N': round(float(geo_info[ii][2]))}).rjust(5)
    string += ('%(ch)0.6f' % {'ch': float(geo_info[ii][1])}).rjust(12)
    for jj in range(3):
      string += ('%(r)0.6f' % {'r': geo_spec[ii][jj]}).rjust(12)
    string += '\n'
  # If exist, write labels
  if labels is not None:
    string += ('%(N)d' % {'N': len(data)}).rjust(5)
    c = 0
    for i,j in enumerate(labels):
      c += 1
      string += str(j).rjust(5)
      if (c  % 9 == 8):
        string += '\n'
    string += '\n'
  # Write data
  for rr in range(len(grid.x)):
    for ss in range(len(grid.y)):
      c = 0
      for tt in range(len(grid.z)):
        for dd in data[:,rr,ss,tt]:
          string += ('%(data).5E' % {'data': dd}).rjust(13)
          if (c % 6 == 5): 
            string += '\n'
          c += 1
      string += '\n'
  
  if filename.endswith('gz'):
    with gzip.open(filename, 'wb') as f:
      f.write(string.encode('utf-8'))
  else:
    with open(filename, 'w') as f:
      f.write(string) def obj_creator(data,filename,geo_info,geo_spec,iso=(-0.01,0.01),**kwargs):
  '''Creates a plain text 3D surface file (obj). 
  
  **Parameters:**
  
  data : numpy.ndarray, shape=N
    Contains the output data.
  filename : str
    Contains the base name of the output file.
  geo_info, geo_spec : 
    See :ref:`Central Variables` for details.
  iso : tuple
    Specifies the isovalue for negative and posititve surfaces, respectively.
  '''
  # Shape shall be (Ndrv,Ndata,Nx,Ny,Nz) or (Ndrv,Ndata,Nxyz)
  data = numpy.array(data)
  dims = 3
  ndim = data.ndim
  shape = data.shape
    
  if data.ndim < dims:
    raise AssertionError('data.ndim < ndim of grid')
  if data.ndim > dims + 1:
    raise AssertionError('data.ndim > (ndim of grid) +2')

  if shape != tuple(grid.N_):
    raise AssertionError('The grid does not fit the data.')

  # Conversion of 3d array and storage as surface file
  basisname = filename.replace('.obj.gz', '').replace('.obj', '')
  spacing = tuple(numpy.concatenate(grid.delta_))
  for i in iso:
        
    # Data conversion    
    verts, faces, normals, values = measure.marching_cubes(data,
                                                           level=i,
                                                           spacing=spacing,
                                                           step_size=1)
    faces = faces + 1
    
    # Define filename for plus and minus surfaces
    pm = 'plus' if i > 0 else 'minus'
    fid = f'{basisname}_{pm}.obj'

    # Write vertices
    string = '# Number of vertices: %d\n' % len(verts)
    for item in verts:
      string += 'v %0.6g %0.6g %0.6g\n' % (item[0],item[1],item[2])
          
    # Write normals
    string += '# Number of normals: %d\n' % len(normals)
    for item in normals:
      string += ('vn %.3g' % item[0]).replace('0.', '.')
      string += (' %0.3g' % item[1]).replace('0.', '.')
      string += (' %0.3g\n' % item[2]).replace('0.', '.')
      
    # Write faces
    string += '# Number of faces: %d\n' % len(faces)
    for item in faces:
      string += 'f {0}//{0} {1}//{1} {2}//{2}\n'.format(item[0],item[1],item[2])

    if filename.endswith('gz'):
      with gzip.open(fid + '.gz', 'wb') as f:
        f.write(string.encode('utf-8'))
    else:
      with open(fid, 'w') as f:
        f.write(string) def hdf5_creator(data, filename, qcinfo=None, gname='', ftype='hdf5', mode='w',attrs={},**kwargs):
  '''Creates an hdf5 file (Hierarchical Data Format) output.
  
  **Parameters:**
  
  data : numpy.ndarray
    Contains the output data.
  filename : str
    Contains the base name of the output file.
  qcinfo : 
    Stores all information from qcinfo class (see :ref:`Central Variables` for
    details).
  gname : str, optional
    Specifies the path, where the data is stored.
  mode : str, optional
    Specifies the mode used to open the file. ('r', 'w', 'a')
  attrs : dict or dict of dicts, optional
    Supply datasets with attributes. If dict, entries will be added to 'data'.
  grid : module or class, global
    Contains the grid, i.e., grid.x, grid.y, and grid.z.
  '''
  
  if ftype.lower() in ['hdf5', 'h5']:
    write = hdf5_write
    attributes = hdf5_attributes
  elif ftype.lower() in ['numpy', 'npz']:
    write = npz_write
    attributes = npz_attributes
  else:
    raise NotImplementedError('File format {0} not implemented for writing.'.format(ftype.lower()))
  
  # Save data and all other kwargs
  write(filename,mode=mode,gname=gname,data=data,**kwargs)
  
  # Save grid
  write(filename,mode='a',gname=os.path.join(gname,'grid'),
             x=grid.x,y=grid.y,z=grid.z,
             is_vector=grid.is_vector,is_regular=grid.is_regular)
  
  # Set attributes
  attributes(filename,gname=gname,**attrs)
  
  if qcinfo is not None:
    write_native(qcinfo, outputname=filename, ftype=ftype, mode='a', 
                 gname=os.path.join(gname,'qcinfo'))   def hdf5_append(x,group,name='data'):
  '''Automatically append data to an open hdf5 file/group.
  
  **Parameters:**
  
  group : h5py.File or h5py.Group
    The hdf5 file/group where the data will be appended.
  name : string, optional
    Specifies the group/dataset name in the hdf5 file/group. If empty,
    root directory of hdf5 file/group is chosen.  
  '''
  from orbkit.orbitals import MOClass, AOClass

  if isinstance(x, MOClass) or isinstance(x, AOClass):
    x = x.todict()
  if isinstance(x,numpy.ndarray):
    if x.dtype.type is numpy.unicode_:
      x = numpy.asarray(x,dtype=numpy.string_)
    h5_dset = group.create_dataset(name,numpy.shape(x),data=x)
  elif isinstance(x,list):
    if name != '':
      subgroup = group.create_group(name)
    else:
      subgroup = group
    num = len(x)
    subgroup.attrs['num'] = num
    keys = range(num)
    for ii in keys:
      hdf5_append(x[ii],subgroup,name='%d' % ii)
  elif isinstance(x,dict):
    if name != '':
      subgroup = group.create_group(name)
    else:
      subgroup = group
    keys = x.keys()
    for ii in keys:
      hdf5_append(x[ii],subgroup,name=ii)
  else:
    if x is None:
      x = str('None')
    group.attrs[name] = x def hdf5_write(fid,mode='w',gname='',**kwargs):
  '''Writes all keyword arguments to an hdf5 file.
  
  **Parameters:**
  
  fid : str
    Contains the filename of the hdf5 file.
  mode : str, optional
    Specifies the mode used to open the file. ('r', 'w', 'a')
  gname : str
    Specifies the group name in the hdf5 file, where the kwargs will be appended
    to.
  
  **Note:**
  
    All additional keyword arguments specify the input data.
  '''
  for f in hdf5_open(fid,mode=mode):
    if gname != '':
      group = f.create_group(gname) if gname not in f.keys() else f[gname]
    else:
      group = f
    for key,data in kwargs.items():
      if isinstance(data,(list,numpy.ndarray)):
        data = numpy.array(data)
        if data.dtype.type is numpy.unicode_:
          data = numpy.asarray(data,dtype=numpy.string_)
        elif data.dtype in ['O']:
          data = numpy.array(data, dtype='S10')
        group[key] = data
      elif isinstance(data,dict):
        hdf5_append(data,group,name=key)
      else:
        if data is None or type(data) == bool:
          data = str(data)
        group.attrs[key] = data def view_with_mayavi(x,y,z,data,is_vectorfield=False,geo_spec=None,datalabels=None,
                     iso_min=1e-4,iso_val=0.01,iso_max=10.0,**kwargs):
  ''' Creates an interactive mayavi dialog showing isosurface plots of the input
  data. 
  
  Components adapted from:
  http://stackoverflow.com/a/1830192
  http://docs.enthought.com/mayavi/mayavi/auto/example_mlab_interactive_dialog.html
  
  **Parameters:**
  
  x,y,z : numpy.ndarray, 1-dim
    Contains the grid.
  data : numpy.ndarray, shape=(len(x), len(y),len(z)) or shape=(N, len(x), len(y),len(z))
    Contains the output data.
  geo_spec : 
    See :ref:`Central Variables` for details.
    If not None, the atom positions will be drawn additionally.
  datalabels : None or list of str
    Contains information about the plotted data with len(datalabels) == len(data).
  '''
  try:
    from enthought.traits.api import HasTraits, Range, Instance, on_trait_change, Bool, Str,List,Button
    from enthought.traitsui.api import View, Item, Group, ListStrEditor,HSplit,VSplit
  except ImportError:
    from traits.api import HasTraits, Range, Instance, on_trait_change, Bool, Str,List,Button
    from traitsui.api import View, Item, Group, ListStrEditor,HSplit,VSplit
  
  try:
    from enthought.mayavi import mlab
    from enthought.mayavi.core.api import PipelineBase
    from enthought.mayavi.core.ui.api import MayaviScene, SceneEditor, MlabSceneModel
  except ImportError:
    from mayavi import mlab
    from mayavi.core.api import PipelineBase
    from mayavi.core.ui.api import MayaviScene, SceneEditor, MlabSceneModel
  
  from copy import deepcopy
  data = numpy.array(data)
  
  if data.ndim == 3:
    data = data[numpy.newaxis]
  elif is_vectorfield and data.ndim == 4:
    data = data[numpy.newaxis]

  if geo_spec is not None and numpy.ndim(geo_spec) == 2:
      geo_spec = [geo_spec for i in range(len(data))]
  
  if datalabels is not None and len(datalabels) < len(data):
    raise ValueError('`datalabels` has to be a list of strings with the same' +
                     'length as `data`.')
  if datalabels is not None:
    datalabels = ['%03d: %s' % (i,j) for i,j in enumerate(datalabels[:len(data)])]
  
  Z,Y,X = meshgrid2(z,y,x)
  
  class MyModel(HasTraits):  
      select  = Range(0, len(data)-1, 0)
      last_select = deepcopy(select)
      iso_value = Range(iso_min, iso_max, iso_val,mode='logslider',
                        label='Iso value' if not is_vectorfield else 'Scale factor')
      if not is_vectorfield:
        opacity = Range(0, 1.0, 1.0)
      else:
        opacity = Range(iso_min, 10*iso_max, 10*iso_val,mode='logslider',
                        label='Range to ')
      show_atoms = Bool(True)
      label = Str()
      available = List(Str)
      available = datalabels
      
      prev_button = Button('Previous')
      next_button = Button('Next')
      
      scene = Instance(MlabSceneModel, ())
      plot_atoms = Instance(PipelineBase)
      plot0 = Instance(PipelineBase)
      
      # When the scene is activated, or when the parameters are changed, we
      # update the plot.
      @on_trait_change('select,iso_value,show_atoms,opacity,label,scene.activated')
      def update_plot(self):        
        if self.plot0 is None:    
          if not is_vectorfield:
            src = mlab.pipeline.scalar_field(X,Y,Z,data[self.select])
            self.plot0 = self.scene.mlab.pipeline.iso_surface(
                        src, contours= [-self.iso_value,self.iso_value], 
                        opacity=self.opacity,colormap='blue-red',
                        vmin=-1e-8,vmax=1e-8)
          else:
            self.plot0 = self.scene.mlab.quiver3d(X,Y,Z,*data[self.select],
                                                  scale_factor=self.iso_value) #flow
          self.plot0.scene.background = (1,1,1)
        elif self.select != self.last_select:
          if not is_vectorfield:
            self.plot0.mlab_source.set(scalars=data[self.select])
          else:
            self.plot0.mlab_source.set(vectors=data[self.select].reshape((3,-1)).T,
                                                  scale_factor=self.iso_value)
          if geo_spec is not None: 
            self.plot_atoms.mlab_source.set(
                                      x=geo_spec[self.select][:,0],
                                      y=geo_spec[self.select][:,1],
                                      z=geo_spec[self.select][:,2])
        if not is_vectorfield:
          self.plot0.contour.contours = [-self.iso_value,self.iso_value]
          self.plot0.actor.property.opacity = self.opacity  
        else:
          #self.plot0.glyph.glyph.range = numpy.array([ 0., self.opacity])
          self.plot0.glyph.glyph.scale_factor = self.iso_value
            
        self.last_select = deepcopy(self.select)
        if datalabels is not None:
          self.label = datalabels[self.select]
        if geo_spec is not None: 
          if self.plot_atoms is None:
            self.plot_atoms = self.scene.mlab.points3d(
                                      geo_spec[self.select][:,0],
                                      geo_spec[self.select][:,1],
                                      geo_spec[self.select][:,2],
                                      scale_factor=0.75,resolution=20)
          self.plot_atoms.visible = self.show_atoms
      
      def _prev_button_fired(self):
        if self.select > 0:
          self.select -= 1
      def _next_button_fired(self):
        if self.select < len(data)-1:
          self.select += 1
      
      
      # The layout of the dialog created
      items = (Item('scene', editor=SceneEditor(scene_class=MayaviScene),
                      height=400, width=600, 
                      show_label=False),)
      items0 = ()
      if len(data) > 1:
        items0 += (Group('select',
                        HSplit(Item('prev_button', show_label=False),
                               Item('next_button', show_label=False)
                              )),)
      items0 += (Group('iso_value', 'opacity','show_atoms'
                        ),)
              
      if datalabels is not None:
        if len(datalabels) > 1:        
            items1 = (Item('available', 
                       editor=ListStrEditor(title='Available Data',editable=False),
                       show_label=False,style='readonly',width=300
                       ),)
            items0 = HSplit(items0,items1)
        items += (Group(Item('label',label='Selected Data',style='readonly', show_label=True),'_'),
               items0,)
      else:
        items += items0
      view = View(VSplit(items[0],items[1:]),
                  resizable=True
                  )
  
  my_model = MyModel()
  my_model.configure_traits() def pdb_creator(geo_info,geo_spec,filename='new',charges=None,comments='',
    angstrom=True):
  '''Creates a plain text file with information concerning the molecular 
  structure in Protein Data Bank (PDB) format. 
  
  **Parameters:**
  
  geo_info, geo_spec : 
    See :ref:`Central Variables` for details.
  filename : str
    Contains the base name of the output file.
  charges : numpy.ndarray, shape=(Natoms,), optional
    Contains a partial charge for each atom.
  comments : str, optional
    Specifies the second (comment) line of the pdb file.
  angstrom : bool, optional
    If True, conversion of molecular coordinates from Bohr radii to Angstrom.
  '''
  
  # Check for charges
  if charges is None:
    charges = numpy.zeros(len(geo_spec))
    
  # Open an empty file 
  fid = open('%(f)s.pdb' % {'f': filename},'w')
  
  # Write HEADER, TITLE and AUTHOR
  fid.write('HEADER\n')
  fid.write('TITLE    %s\n' % comments)
  fid.write('AUTHOR    orbkit\n') 
  
  # Write ATOM records
  string = ''
  for il in range(len(geo_spec)):
    string = 'ATOM    %s' % (il+1)
    while len(string) < 13:
      string += ' '
    string += '%s' % (geo_info[il][0])
    while len(string) < 17:
      string += ' '
    xyz = geo_spec[il]
    if angstrom:
      xyz *= a0_to_aa
    string += '             %s %s %s        ' % (('%.9f' % xyz[0])[:7],
                                                 ('%.9f' % xyz[1])[:7],
                                                 ('%.9f' % xyz[2])[:7])
    string += '%s        ' % (('%.9f' % charges[il])[:6])
    fid.write(str(string) + '\n')
      
  # Write MASTER and END line
  fid.write('MASTER        0    0    0    0    0    0    0    0 ' + 
            '%s    0    0    0\nEND' % ('%s'.rjust(4) % len(geo_spec)))
  
  # Close the file 
  fid.close() def vmd_network_creator(filename,cube_files=None,render=False,exit=False,
                        iso=(-0.01,0.01),abspath=False,mo_options='',**kwargs):
  '''Creates a VMD script file from a list of cube files provided.

  **Parameters:**

  filename : str
    Contains the base name of the output file.
  cube_files : None or list of str
    Specifies the cube files which serve as input for the VMD script.
    If None, searches the directory for '.cb' and '.cube' files.
  render : bool
    If True, the VMD script will automatically create '.tga' files for each
    cube file.
  exit : bool
    If True, exit vmd after rendering.
  iso : tuple
    Specifies the isovalue for the blue and the red isosurface, respectively.
  abspath : bool
    If True, the paths of the cube files will be expanded to absolute file paths.
  mo_options : strings
    Further commands to include in vmd script.
    
  **Example:**
  
  .. code-block::
  
    # create cube files
    extras.calc_mo(qc, active, otype=['vmd'], ofid=fname)

    # recreate vmd script
    mo_options = """
    rotate y by 90
    rotate x by 45
    rotate z by 0
    axes location off
    display distance 1
    """
    
    vmd_network_creator(fname, render=1, exit=1, mo_options=mo_options)
  '''
  from os import path,listdir
  import linecache
  from . import vmd_network_draft
  
  filename += '.vmd' if not filename.endswith('vmd') else ''
  
  if cube_files is None:
    display('No list of cube (.cb or .cube) filenames provided. Checking the directory' +
            ' of the outputfile...')
    cube_files = []
    for fid in listdir(path.dirname(filename)):
      if fid.endswith('.cb') or fid.endswith('.cube'):
        cube_files.append(fid)
    if cube_files == []:
      raise IOError('Could not find valid cube files in %s' % path.dirname(filename))
  elif isinstance(cube_files,str):
    cube_files = [cube_files]
  elif not isinstance(cube_files,list):
    raise IOError('`cube_files` has to be a list of strings.')

  title = []
  mo = ''
  for i,f in enumerate(cube_files):
    title = linecache.getline(f,2)
    if title.split() == []:
      title = path.splitext(path.basename(f))[0]
    else:
      title = title.replace('\n','').replace(' ','')
    linecache.clearcache()
    pid = path.abspath(f) if abspath else path.relpath(f,path.dirname(filename))
    mo += vmd_network_draft.mo_string % {
                                  'c': i,
                                  'n1': pid,
                                  'n2': title,
                                  'isored': iso[0],
                                  'isoblue': iso[1],
                                  'render': '' if render else '#',
                                  'mo_options' : mo_options,
                                  }
    
  with open(filename,'w') as f:
    f.write(vmd_network_draft.vmd_string % {'mo':mo})
    if exit:
      f.write('exit\n')
  f.close() def xyz_creator(geo_info,geo_spec,filename='new',mode='w',charges=None,comments='',
    angstrom=True):
  '''Creates a xyz file containing the molecular coordinates. 
  
  **Parameters:**
  
  geo_info, geo_spec : 
    See :ref:`Central Variables` for details.
  filename : str
    Contains the base name of the output file.
  charges : numpy.ndarray, shape=(Natoms,), optional
    Contains a partial charge for each atom.
  comments : str, optional
    Specifies the second (comment) line of the xyz file.
  angstrom : bool, optional
    If True, conversion of molecular coordinates from Bohr radii to Angstrom.
  '''
  
  # Open an empty file 
  fid = open('%s.xyz' % filename,mode)
  
  # Write number of atoms and a comment line
  fid.write('%d\n%s\n' % (len(geo_spec),comments))
  
  # Write Cartesian coordinates of molecular structure
  string = ''
  for il in range(len(geo_spec)):
    string += '%-2s' % geo_info[il][0]
    for i in range(3):
      xyz = geo_spec[il]
      if angstrom:
        xyz *= a0_to_aa
      string += ' %22.15f'  % (xyz[i])
    if charges is not None:
      string += ' %22.15f'  % charges[il]
    string += '\n'
  fid.write(string)
  
  # Close the file
  fid.close() def write_native(outdata, outputname, ftype='numpy', gname='qcinfo', mode='w',
                 **add_data):
  '''Creates the requested output.
  
  **Parameters:**

    outdata : 
      Instance of the class to be written to fine. This must support a **todict** function.
      outputname : str or list of str
      Contains the base name of the output file.
    ftype, str : 
      Data can be written to ``.npz`` or ``hdf5`` files.
  '''

  if not callable(getattr(outdata, 'todict')):
    raise NotImplementedError('Supplied class instance does not support todict()')

  if ftype.lower() in ['hdf5', 'h5']:
    from . import hdf5_write
    write = hdf5_write
    if not (outputname.endswith('hdf5') or outputname.endswith('h5')):
      outputname += '.hdf5'
  elif ftype.lower() in ['numpy', 'npz']:
    from . import npz_write
    write = npz_write
    if not (outputname.endswith('numpy') or outputname.endswith('npz')):
      outputname += '.npz'
  else:
    raise NotImplementedError('File format {0} not implemented for writing.'.format(ftype.lower()))
    #raise NotImplementedError('Only .npz and .hdf5 are currently supported.')
  
  odata = outdata.todict()
  odata.update(add_data)
  odata['date'] = time.strftime("%Y-%m-%d") 
  odata['time'] = time.strftime("%H:%M:%S")
  odata = unravel_dicts(odata)
  for i,key in enumerate(odata.keys()):
    if i == 0:
      current_mode = mode
    else:
      current_mode = 'a'
    write(outputname, mode=current_mode, gname=join(gname,key), **odata[key])
  
  if gname != '':
    return '{0}@{1}'.format(outputname,gname)
  else:
    return outputname
import numpy
import h5py
from os import path
from orbkit import grid, options
from orbkit.display import display
from orbkit.units import a0_to_aa
from .amira import amira_creator, hx_network_creator
from .cube import cube_creator
from .obj import obj_creator
from .hdf5 import hdf5_creator, hdf5_append, hdf5_write
from .mayavi_interface import view_with_mayavi
from .pdb import pdb_creator
from .vmd import vmd_network_creator
from .xyz import xyz_creator
from .native import write_native
the License, or any later version.

orbkit 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public 
License along with orbkit.  If not, see <http://www.gnu.org/licenses/>.
'''

# Import general modules

# Import orbkit modules


synonyms = {'auto':'auto',
            'h5':'h5', 'hdf5':'h5',
            'npz':'npz','numpy':'npz',
            'cube':'cube', 'cb':'cube',
            'cube.gz':'cube', 'cb.gz':'cube',
            'obj':'obj', 'obj.gz':'obj',
            'am':'am', 
            'hx':'hx',
            'vmd':'vmd',
            'mayavi':'mayavi',
            '':None
            }

def main_output(data,qc=None,outputname='data',otype='auto',gname='',
                drv=None,omit=[],datalabels='',dataindices=None,mode='w',**kwargs):
  '''Creates the requested output.
  
  **Parameters:**
  
  data : numpy.ndarray, shape=N, shape=((NDRV,) + N), shape=(n, (NDRV,) + N) or list of numpy.ndarrays
    Contains the output data. The shape (N) depends on the grid and the data, i.e.,
    3d for regular grid, 1d for vector grid. 
  qc : class or dict
    QCinfo class or dictionary containing the following attributes/keys.
    See :ref:`Central Variables` for details.
  outputname : str or list of str
    Contains the base name of the output file. If outputname contains @, string will be split and first 
    part interpreted as outputname and second as gname (cf. Parameters:gname). 
  otype : str or list of str, optional
    Contains the output file type. Possible options:
    'auto', 'h5', 'cb', 'am', 'hx', 'vmd', 'mayavi', 'obj'
    
    If otype='native', a native input file will be written, the type of which may be specifie by
    ftype='numpy'.
  gname : str, optional
    For native, HDF5, or npz output, specifies the group, where the data will be stored.
  drv : None, list of str or list of list of str, optional
    If not None, a 4d(regular)/2d(vector) input data array will be expected
    with NDRV = len(drv). Specifies the file labels, i.e. e.g., data_d{drv}.cube for 4d array.
    For 5d arrays i.e., data_0_d{drv}.cube
  datalabels : list of str, optional
    If not empty, the output file types specified here are omitted.
  dataindices : list of int, optional
    If not empty, the output file types specified here are omitted.
  omit : list of str, optional
    If not empty, the output file types specified here are omitted.
  mode : str={'r', 'w', 'a'}, optional
    Specifies the mode used to open the file (native, HDF5, or npz). 
  
  **Note:**
  
    All additional keyword arguments are forwarded to the output functions.
  '''
  if otype is None or otype == []:
    return []
  
  if isinstance(outputname, str):
    if '@' in outputname:
      outputname,gname = outputname.split('@')
  if isinstance(otype, str):
    if otype == 'auto':
      outputname, otype = path.splitext(outputname)
      otype = otype[1:]
    otype = [otype]
  elif isinstance(otype, list) and len(otype) == 1:
    if otype[0] == 'auto':
      outputname, otype = path.splitext(outputname)
      otype = [otype[1:]]
  else:
    for iot in range(len(otype)):
      if otype[iot] == 'auto':
        outputname, tmp = path.splitext(outputname)
        if tmp != '':
          otype[iot] = tmp[1:]
  
  # Catch our native format before all else
  # We can't figure this out by the file ending alone
  # as we support hdf5 for output of both grid-based data
  # as well as our internal format

  output_written = []
  internals = [i for i in range(len(otype)) if otype[i] == 'native']
  
  if len(internals) > 0:
    if not isinstance(data,list):
      data = [data]

    if isinstance(outputname,str):
      outputname = [outputname for _ in data]

    if 'ftype' in kwargs.keys():
      if isinstance(kwargs['ftype'],str):
        ftype = [kwargs['ftype'] for _ in data]
    else:
      ftype = ['numpy' for _ in data]

    if 'group' in kwargs.keys():
      if isinstance(kwargs['group'],str):
        group = [kwargs['group'] for _ in data]
    else:
      group = [i.__class__.__name__.lower() for i in data]
    
    display('Writing native input file...' )
    for i, oname in enumerate(outputname):

      output_written.append(write_native(data[i], oname, ftype[i], mode=mode,