# Path: orbkit/read/high_level.py
# def main_read(fname, all_mo=False, spin=None, itype='auto', check_norm=False, **kwargs):
#   '''
#   This is the high-lever interface for the
#   orbkit reading routines.
#   
#   **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.
#   all_mo : bool, optional
#     If True, all molecular orbitals are returned.
#   spin : {None, 'alpha', or 'beta'}, optional
#     If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
#   itype : str, optional
#     Can be used to manually specify the input filetype.
#   check_norm : bool, optional
#     If True, ORBKIT verifies that molecular orbitals are orthonormal.
#     
#   **Note:**
#   
#     All additional keyword arguments are forwarded to the reading functions.
#     
#   **Returns:**
#   
#     qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
#     See :ref:`Central Variables` for details.
#   '''
# 
#   if isinstance(fname, str):
#     filename = fname
#   else:
#     filename = fname.name
# 
#   if itype == 'auto':
#     itype = find_itype(fname)
#  
#   display('Loading data from {0} type file {1}\n'.format(itype, filename))
# 
#   qc = readers[itype](fname, all_mo=all_mo, spin=spin, **kwargs)
# 
#   if check_norm:
#     deviation = check_mo_norm(qc)
#     if deviation >= 1e-5:
#       raise ValueError('Bad molecular orbital norm: {0:%4e}'.format(deviation))
# 
#   return qc
# 
# Path: orbkit/analytical_integrals.py
# def get_ao_overlap(coord_a, coord_b, ao_spec, lxlylz_b=None,
#                    drv=None):
# def cartesian2spherical_aoom(ao_overlap_matrix,ao_spec):
# def get_mo_overlap(mo_a,mo_b,ao_overlap_matrix):
# def initializer(gargs):
# def get_slice(x):  
# def get_mo_overlap_matrix(mo_a,mo_b,ao_overlap_matrix,numproc=1):
# def get_moom_atoms(atoms,qc,mo_a,mo_b,ao_overlap_matrix,numproc=1):
# def check_mo_norm(qc):
# def get_dipole_moment(qc,component=['x','y','z']):
# def get_ao_dipole_matrix(qc,component='x'):
# def get_nuclear_dipole_moment(qc,component='x'):
# def get_atom2mo(qc):
# def get_lc(atoms,atom2mo,strict=False):


from orbkit.read import main_read
from orbkit.tools import *
from orbkit.analytical_integrals import get_ao_overlap,get_mo_overlap,print2D

print('The total number of electrons is %.8f' % electron_number)
# Compute the x-, y-, and z-component of the dipole moment
dipole_moment = []
for component in range(3):
  
  # Compute the first part of the expectation value:
  # Get the the exponents lx, ly, lz for the primitive Cartesian Gaussians of
  # the `Ket` basis set, and increase lz by one.
  lxlylz_b = qc.ao_spec.get_lxlylz()
  lxlylz_b[:,component] += 1

  ao_part_1 = get_ao_overlap(qc.geo_spec,qc.geo_spec,qc.ao_spec,
                             lxlylz_b=lxlylz_b)

  # Compute the second part of the expectation value:
  ao_part_2 = get_ao_overlap(qc.geo_spec,qc.geo_spec,qc.ao_spec)

  i = 0
  for sel_ao in range(len(qc.ao_spec)):
    l = l_deg(l=qc.ao_spec[sel_ao]['type'].lower(),
              cartesian_basis=(not qc.ao_spec.spherical))
    for ll in range(l):
      ao_part_2[:,i] *= qc.geo_spec[qc.ao_spec[sel_ao]['atom'],component]
      i += 1
  
  ao_dipole_matrix = (ao_part_1+ao_part_2)
  
  # Print the atomic orbital dipole matrix
  if 0:

# Path: orbkit/read/high_level.py
# def main_read(fname, all_mo=False, spin=None, itype='auto', check_norm=False, **kwargs):
#   '''
#   This is the high-lever interface for the
#   orbkit reading routines.
#   
#   **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.
#   all_mo : bool, optional
#     If True, all molecular orbitals are returned.
#   spin : {None, 'alpha', or 'beta'}, optional
#     If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
#   itype : str, optional
#     Can be used to manually specify the input filetype.
#   check_norm : bool, optional
#     If True, ORBKIT verifies that molecular orbitals are orthonormal.
#     
#   **Note:**
#   
#     All additional keyword arguments are forwarded to the reading functions.
#     
#   **Returns:**
#   
#     qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
#     See :ref:`Central Variables` for details.
#   '''
# 
#   if isinstance(fname, str):
#     filename = fname
#   else:
#     filename = fname.name
# 
#   if itype == 'auto':
#     itype = find_itype(fname)
#  
#   display('Loading data from {0} type file {1}\n'.format(itype, filename))
# 
#   qc = readers[itype](fname, all_mo=all_mo, spin=spin, **kwargs)
# 
#   if check_norm:
#     deviation = check_mo_norm(qc)
#     if deviation >= 1e-5:
#       raise ValueError('Bad molecular orbital norm: {0:%4e}'.format(deviation))
# 
#   return qc
# 
# Path: orbkit/analytical_integrals.py
# def get_ao_overlap(coord_a, coord_b, ao_spec, lxlylz_b=None,
#                    drv=None):
# def cartesian2spherical_aoom(ao_overlap_matrix,ao_spec):
# def get_mo_overlap(mo_a,mo_b,ao_overlap_matrix):
# def initializer(gargs):
# def get_slice(x):  
# def get_mo_overlap_matrix(mo_a,mo_b,ao_overlap_matrix,numproc=1):
# def get_moom_atoms(atoms,qc,mo_a,mo_b,ao_overlap_matrix,numproc=1):
# def check_mo_norm(qc):
# def get_dipole_moment(qc,component=['x','y','z']):
# def get_ao_dipole_matrix(qc,component='x'):
# def get_nuclear_dipole_moment(qc,component='x'):
# def get_atom2mo(qc):
# def get_lc(atoms,atom2mo,strict=False):


# Path: examples/basic_examples/orbkit_non-grid-based.py
from orbkit.read import main_read
from orbkit.tools import *
from orbkit.analytical_integrals import get_ao_overlap,get_mo_overlap,print2D


print('The total number of electrons is %.8f' % electron_number)
# Compute the x-, y-, and z-component of the dipole moment
dipole_moment = []
for component in range(3):
  
  # Compute the first part of the expectation value:
  # Get the the exponents lx, ly, lz for the primitive Cartesian Gaussians of
  # the `Ket` basis set, and increase lz by one.
  lxlylz_b = qc.ao_spec.get_lxlylz()
  lxlylz_b[:,component] += 1

  ao_part_1 = get_ao_overlap(qc.geo_spec,qc.geo_spec,qc.ao_spec,
                             lxlylz_b=lxlylz_b)

  # Compute the second part of the expectation value:
  ao_part_2 = get_ao_overlap(qc.geo_spec,qc.geo_spec,qc.ao_spec)

  i = 0
  for sel_ao in range(len(qc.ao_spec)):
    l = l_deg(l=qc.ao_spec[sel_ao]['type'].lower(),
              cartesian_basis=(not qc.ao_spec.spherical))
    for ll in range(l):
      ao_part_2[:,i] *= qc.geo_spec[qc.ao_spec[sel_ao]['atom'],component]
      i += 1
  
  ao_dipole_matrix = (ao_part_1+ao_part_2)
  
  # Print the atomic orbital dipole matrix
  if 0:

    print2D(ao_dipole_matrix) 