def alignSVD(self, mol, ref_list=None, tar_list=None):
    if type(mol) is str:
      try:
        mol = qtk.Molecule(mol)
      except:
        qtk.exit("error when reading molecule file: %s" % mol)
    assert issubclass(mol.__class__, qtk.Molecule)
    if not ref_list:
      ref_list = [i for i in range(self.N)]
    if not tar_list:
      tar_list = copy.deepcopy(ref_list)

    lst_a = self.R[ref_list]
    lst_b = mol.R[tar_list]
    center_a = np.mean(lst_a, axis=0)
    center_b = np.mean(lst_b, axis=0)
    #na = len(lst_a)
    na = self.N
    #nb = len(lst_b)
    nb = mol.N
    crd_a = self.R - np.kron(center_a, np.ones((self.N, 1)))
    crd_b = mol.R - np.kron(center_b, np.ones((mol.N, 1)))
    ref_a = lst_a - np.kron(center_a, np.ones((len(lst_a), 1)))
    ref_b = lst_b - np.kron(center_a, np.ones((len(lst_b), 1)))

    H = np.dot(np.transpose(ref_a), ref_b)
    U, s, V = np.linalg.svd(H)
    R = np.dot(np.transpose(V), np.transpose(U))
    self.R = np.transpose(
      np.dot(R, np.transpose(crd_a))) + \
      np.kron(center_b, np.ones((na, 1))
    )

  def sort(self, order = 'Zxyz'):
    odict = {'x':0, 'y':1, 'z':2}
    tmp = []
    for o in order:
      if o == 'Z':
        tmp.insert(0, self.Z)
      elif o in odict:
        tmp.insert(0, self.R[:, odict[o]])
      else:
        qtk.exit("sorting order '%c' not valid" % o)
    ind = np.lexsort(tmp)
    self.R = self.R[ind]
    if list(self.R_scale[0]):
      self.R_scale = self.R_scale[ind]
    self.Z = list(np.array(self.Z)[ind])
    self.type_list = list(np.array(self.type_list)[ind])
    self.string = list(np.array(self.string)[ind])

    if order == 'Zxyz':
      type_list = []
      self.index = []
      for i in range(self.N):
        Zn = self.type_list[i]
        if Zn not in type_list:
          type_list.append(Zn)
          self.index.append(i)
      self.index.append(self.N)

 def setReference(self, ref_vec):
   if len(ref_vec) != self.data_size:
     qtk.exit("number of data points not match")
   def set_ref(i, ref):
     self.data[i].ref = ref
   vset_ref = np.vectorize(set_ref)
   vset_ref(range(self.data_size), ref_vec)

  def __init__(self, *args, **kwargs):
    if len(args) == 1 and args[0] > 0:
      self.base_dim = args[0]
    elif len(args) == 0 or args[0] == 0:
      fxyz = open(kwargs['xyz'], "r")
      self.base_dim = int(fxyz.readline())
      fxyz.close()
    
    if 'xyz' in kwargs:
      self.data = cm.coulomb_matrix(kwargs['xyz'], 
                                         self.base_dim)
      self.name = re.sub('.*\/','',kwargs['xyz'])
      self.name = re.sub('\.xyz','',self.name)
    else:
      qtk.exit("CoulombMatrix: input mode is not specified")

    # set row norm matrix
    NORM = np.vstack([sum(self.data), range(self.base_dim)]).T
    # sort NORM by row norm
    NORM=NORM[NORM[:,0].argsort()]
    # reverse array order
    NORM=NORM[::-1]
    # extract new row order
    sortRow=NORM[:,1].astype(int)
    # rearrange Coulomb matrix
    self.data = self.data[:,sortRow][sortRow]

  def brillouinize(self, kmesh=None, kgrid=None):

    if not hasattr(self, 'kpoints') or len(self.kpoints) == 0:
      qtk.exit('kpoints information not available')
    else:
      if not hasattr(self, 'kpoints_symmetrized'):
        k_old = self.kpoints[:,:3]
        b_old = self.band
      else:
        k_old = self.kpoints_symmetrized[:,:3]
        b_old = self.band_symmetrized

    if kgrid is None and kmesh is not None:
      try:
        kgrid, new_band = self._spg_grid(k_old, b_old, kmesh)
      except:
        try:
          kgrid, new_band = self._ase_grid(k_old, b_old, kmesh)
        except Exception as e:
          qtk.exit('kpoint generation failed')
    elif kmesh is None and kgrid is not None:
      new_band = self._kgrid_template(kgrid)

    if not hasattr(self, 'kpoints_symmetrized'):
      self.kpoints_symmetrized = self.kpoints
      self.band_symmetrized = self.band
    self.kpoints = kgrid
    self.band = new_band

def coulomb_matrix(mol, n = -1, size = 0, 
                   sort = True, nuclear_charges = True):
  if size == 0:
    size = mol.N
  if size < mol.N:
    qtk.exit("matrix size too small")
  positions = mol.R
  if nuclear_charges:
    charges = np.array(mol.Z)
  else:
    charges = np.ones(mol.N)
  differences = positions[:, np.newaxis, :] \
              - positions[np.newaxis, :, :]
  distances = np.sqrt((differences ** 2).sum(axis=-1))
  distances[distances == 0] = np.nan # replace 0 for division
  if n != 0:
    invR = (distances ** n)
  else:
    invR = distances
  invR[np.isnan(invR)] = 0 # change 0 back for getting diagonal
  diag_mask = (invR == 0).astype(int)
  charge_mask_Zij = charges[:, np.newaxis] \
                  * charges[np.newaxis, :]
  charge_mask_2p4 = 0.5 * ((charges[:, np.newaxis] \
                            * charges[np.newaxis, :]) \
                            * diag_mask) ** 1.2
  cm = invR * charge_mask_Zij + charge_mask_2p4
  if sort:
    ind = np.argsort(cm.sum(axis=-1))
    cm = cm[:, ind][ind]
  out = np.zeros([size, size])
  out[:cm.shape[0], :cm.shape[1]] = cm
  return out

 def _kPath(self, special_points, dk):
   spk = self.special_kpoints
   out = []
   pos = 0
   tick_pos = [0]
   tick_txt = []
   for i in range(len(special_points)-1):
     ci = special_points[i]
     cf = special_points[i+1]
     if ci not in spk:
       qtk.exit("special point %c not in reconized" % ci)
     if cf not in spk:
       qtk.exit("special point %c not in reconized" % cf)
     if len(tick_txt) == 0:
       tick_txt.append("$\mathrm{" + ci + "}$")
     Ri = spk[ci]
     Rf = spk[cf]
     if len(out) == 0:
       out.append(Ri)
     vector = np.array(Rf) - np.array(Ri)
     n_points = int(ceil(np.linalg.norm(vector) / float(dk)))
     pos = pos + n_points
     tick_pos.append(pos)
     tick_txt.append("$\mathrm{" + cf + "}$")
     for i in range(1, n_points+1 ):
       point = list(np.round(
         np.array(Ri) + (i/float(n_points))*vector, decimals=3
       ))
       out.append(point)
     tick_txt = [
       '$\Gamma$' if x == '$\mathrm{G}$' else x for x in tick_txt]
   return out, [tick_pos, tick_txt]

  def mo_g09_nwchem(self):
    mo = self.mo_vectors
    ind = np.arange(len(mo[0]))
    itr = 0

    # hard coded reordering for d and f orbitals
    order = {
     'd': [0, 3, 4, 1, 5, 2],
     'f': [0, 4, 5, 3, 9, 6, 1, 8, 7, 2],
    }

    while itr < len(mo[0]):
      bStr = self.basis[itr]['type']
      if len(bStr) > 2:
        key = bStr[0]
        if key in order.keys():
          ind_itr = ind[itr: itr+len(order[key])]
          ind_itr = ind_itr[order[key]]
          for i in range(len(order[bStr[0]])):
            ind[itr + i] = ind_itr[i]
          itr += len(order[bStr[0]]) - 1
        else:
          qtk.exit("basis reordering for %s orbital " % key\
                 + "not implemented yet")
      itr += 1
    return mo[:, ind]

 def cm_check(self, mol):
   ve = mol.getValenceElectrons()
   if (ve % 2) == (mol.multiplicity % 2):
     msg = "Multiplicity %d " % mol.multiplicity + \
           "and %d valence electrons " % ve +\
           "\n(with charge %3.1f) " % float(mol.charge) +\
           "are not compatible"
     qtk.exit(msg)

def n2m(Zn):
    match = [m for m in mass_list.iterkeys() if m in Zn]
    mlen = [len(s) for s in match]
    if len(match) > 0:
        ind = np.argmax(mlen)
        return float(mass_list[match[ind]])
    else:
        qtk.exit("n2Z: element type " + str(Zn) + " is not defined")

def convE(source, units, separator=None):
  def returnError(ioStr, unitStr):
    msg = 'supported units are:\n'
    for key in Eh.iterkeys():
      msg = msg + key + '\n'
    qtk.report(msg, color=None)
    qtk.exit(ioStr + " unit: " + unitStr + " is not reconized")

  EhKey = {
    'ha': 'Eh',
    'eh': 'Eh',
    'hartree': 'Eh',
    'ry': 'Ry',
    'j': 'J',
    'joule': 'J',
    'kj/mol': 'kJ/mol',
    'kjmol': 'kJ/mol',
    'kjm': 'kJ/mol',
    'kj': 'kJ/mol', # assume no kilo Joule!
    'kcal/mol': 'kcal/mol',
    'kcalmol': 'kcal/mol',
    'kcal': 'kcal/mol', # assume no kcal!
    'kcm': 'kcal/mol',
    'ev': 'eV',
    'cminv': 'cmInv',
    'cminverse': 'cmInv',
    'icm': 'cmInv',
    'cm-1': 'cmInv',
    'k': 'K',
    'kelvin': 'K',
  }
 
  Eh = {
    'Eh': 1.0,
    'Ry': 2.0,
    'eV': 27.211396132,
    'kcal/mol': 627.509469,
    'cmInv': 219474.6313705,
    'K': 3.15774646E5,
    'J': 4.3597443419E-18,
    'kJ/mol': 2625.49962
  }

  if not separator:
    separator='-'
  unit = units.split(separator)
  if len(unit) != 2:
    qtk.exit("problem with unit separator '%s'" % separator)
  if unit[0].lower() != 'hartree' and unit[0].lower() != 'eh':
    if unit[0].lower() in EhKey:
      unit0 = EhKey[unit[0].lower()]
      source = source / Eh[unit0]
    else: returnError('input', unit[0])
  if unit[1].lower() not in EhKey: 
    returnError('output', unit[1])
  else:
    unit1 = EhKey[unit[1].lower()]
  return source * Eh[unit1], unit1

def qAtomicNumber(query):
    if type(query) == str:
        if z_list.has_key(query):
            return n2Z(query)
    elif type(query) == int or type(query) == float:
        if type_list.has_key(int(query)):
            return query
    else:
        qtk.exit("qAtom: element " + str(Zn) + " is not defined")

  def __init__(self, molecule, **kwargs):

    if not found:
      qtk.exit("horton module not found.")
    if 'wf_convergence' not in kwargs:
      kwargs['wf_convergence'] = 1e-06
    GaussianBasisInput.__init__(self, molecule, **kwargs)
    self.setting.update(kwargs)
    self.backup()

    mol = IOData(coordinates=molecule.R, numbers=molecule.Z)
    obasis = get_gobasis(mol.coordinates, mol.numbers,
                         self.setting['basis_set'])
    grid = BeckeMolGrid(mol.coordinates, mol.numbers, 
                        mol.pseudo_numbers)

    # Create a linalg factory
    lf = DenseLinalgFactory(obasis.nbasis)
    
    # Compute Gaussian integrals
    olp = obasis.compute_overlap(lf)
    kin = obasis.compute_kinetic(lf)
    na = obasis.compute_nuclear_attraction(mol.coordinates, 
                                           mol.pseudo_numbers, lf)
    er = obasis.compute_electron_repulsion(lf)

    # Create alpha orbitals
    exp_alpha = lf.create_expansion()
    
    # Initial guess
    guess_core_hamiltonian(olp, kin, na, exp_alpha)

    external = {'nn': compute_nucnuc(mol.coordinates, 
                                     mol.pseudo_numbers)}
    libxc_term = RLibXCHybridGGA('xc_b3lyp')
    terms = [
        #RTwoIndexTerm(kin, 'kin'),
        RDirectTerm(er, 'hartree'),
        RGridGroup(obasis, grid, [libxc_term]),
        RExchangeTerm(er, 'x_hf', libxc_term.get_exx_fraction()),
        RTwoIndexTerm(na, 'ne'),
    ]
    ham = REffHam(terms, external)

    self.ht_mol = mol
    self.ht_grid = grid
    self.ht_external = external
    self.ht_obasis = obasis
    self.ht_lf = lf
    self.ht_olp = olp
    self.ht_kin = kin
    self.ht_na = na
    self.ht_er = er
    self.ht_exp_alpha = exp_alpha
    self.ht_ham = ham

def Z2n(Z):
    if type_list.has_key(Z):
        return type_list[Z]
    elif type(Z) is float:
        return "ATOM_%4.2f" % Z
        qtk.warning("Z2n: atomic number not defined, return HETATM")
    elif type(Z) is int:
        return "ATOM_%d" % Z
        qtk.warning("Z2n: atomic number not defined, return HETATM")
    else:
        qtk.exit("Z2n: atomic number " + str(Z) + " is not defined")

  def __init__(self, penalty, penalty_input, inpgen, 
               mating_function, pop_size, **kwargs):
    # redefine some functions are necessary!
    opt.Optimizer.__init__(self, penalty, penalty_input, 
                           inpgen, **kwargs)
    self.pop_size = pop_size
    self.mating_function = mating_function
    self.mutation_rate = 0.05

    if not pop_size/2 > 2:
      print pop_size
      qtk.exit("population too small")

 def loadCubeList(cls, path_list, program=qtk.setting.qmcode):
   if program == 'cpmd':
     cls._file_list = path_list
     _para = [[name] for name in cls._file_list]
     if len(_para)<3*qtk.setting.cpu_count:
       cls._cube_list = qtk.parallelize(qtk.CUBE, _para, 
                                    block_size=1)
     else:
       cls._cube_list = qtk.parallelize(qtk.CUBE, _para)
   else:
     qtk.exit("density of alchemical path is "\
              +"not yet implemented for %s" % self.program)

  def __init__(self, *args, **kwargs):
    # number of atoms
    self.N = 0
    # atom coordinates
    self.R = np.atleast_2d(np.array([]))
    self.R_scale = np.atleast_2d(np.array([]))
    # atom symbols
    self.type_list = []
    # nuclear charges
    self.Z = []
    # moelcule charge
    self.charge = 0
    self.multiplicity = 1
    # index of different atoms
    self.index = 0
    self.bonds = {}
    self.bond_types = {}
    self.string = []
    self.segments = []
    self.periodic = False
    self.scale = False
    self.celldm = False
    self.symmetry = False
    self.grid = False
    self.name = ''

    if len(args) == 1:
      mol = args[0]
      self.read(mol, **kwargs)
    elif len(args) == 2:
      N = len(args[0])
      dim1 = np.array(args[0]).shape
      dim2 = np.array(args[1]).shape
      if dim1 == (N,) and dim2 == (N, 3):
        atoms = args[0]
        coord = np.array(args[1])
      elif dim1 == (N, 3) and dim2 == (N,):
        atoms = args[1]
        coord = np.array(args[0])
      else:
        qtk.exit('not supported declaration of molecule object.')
      self.addAtoms(atoms, coord)
     

    attr_list = dir(self)
    for string, value in kwargs.iteritems():
      if string in attr_list:
        setattr(self, string, kwargs[string])

    if self.N > 0:
      self.ve = self.getValenceElectrons()
      self.ne = sum(self.Z)

 def getRho(self, cartesian=True, resolution='fine', new=False, occupation=None, **kwargs):
   if 'gridpoints' in kwargs:
     new = True
   if new or not hasattr(self, '_rho'):
     self.getPsi(cartesian, resolution, new, **kwargs)
     if not hasattr(self, 'occupation'):
       qtk.exit("occupation number not found")
     if occupation is None:
       occ = np.array(self.occupation)
     else:
       occ = np.array(occupation)
     self._rho = np.sum(self._psi**2 * occ[:, np.newaxis], axis = 0)
   return self._rho

 def __add__(self, other):
   if self.R_scale.shape[1] > 0:
     qtk.exit('Molecule add not implemented for crystals.' + \
              'use extend/setAtoms instead')
   out = Molecule()
   out.N = self.N + other.N
   out.R = np.vstack([self.R, other.R])
   out.Z = np.hstack([self.Z, other.Z])
   out.type_list = np.hstack([self.type_list, other.type_list])
   out.string = np.hstack([self.string, other.string])
   out.charge = self.charge + other.charge
   out.name = self.name + "_" + other.name
   return out

 def getDPsi(self, cartesian=True, resolution='fine', new=False, **kwargs):
   if 'gridpoints' in kwargs:
     new = True
   if new or not hasattr(self, '_dpsi'):
     self.getDPhi(cartesian, resolution, new, **kwargs)
     if not hasattr(self, 'mo_vectors'):
       qtk.exit('mo_vectors not found')
     mo = self.mo_vectors
     if hasattr(self, 'program'):
       if self.program == 'gaussian':
         mo = self.mo_g09_nwchem()
     self._dpsi = np.dot(mo, np.swapaxes(self._dphi, 0, 1))
   return self._dpsi