# 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

First Part:

  Compute the electron number analytically using the standard density formula.

Second Part:

  Compute the dipole moments analytically. This expectation value is more complex.
  As an explanation, consider the following expectation value between two primitive 
  $p_z$ orbitals:

  .. math::

    <\phi_k|z|\phi_l> = \int dxdydz (z - Z_k)^{k_z} e^{\alpha_k r_k^2} (z - Z_l)^{l_z+1} e^{\alpha_l r_l^2} 
                        + Z_l \int dxdydz (z - Z_k)^{k_z} e^{\alpha_k r_k^2} (z - Z_l)^{l_z} e^{\alpha_l r_l^2} 
                      = ao_part_1 + ao_part_2

with :math:`r_l = \sqrt{(x-X_l)^2 + (y-Y_l)^2 + (z - Z_l)^2)}`
'''

in_fid = 'h2o.molden'
# Read the input file
qc = main_read(in_fid,itype='molden',all_mo=False)

# Compute atomic orbital overlap matrix 
ao_overlap_matrix = get_ao_overlap(qc.geo_spec,qc.geo_spec,qc.ao_spec)

# Compute the overlap of the molecular orbitals and weight it with the occupation number
electron_number = 0.
for i_mo in qc.mo_spec:

# 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


First Part:

  Compute the electron number analytically using the standard density formula.

Second Part:

  Compute the dipole moments analytically. This expectation value is more complex.
  As an explanation, consider the following expectation value between two primitive 
  $p_z$ orbitals:

  .. math::

    <\phi_k|z|\phi_l> = \int dxdydz (z - Z_k)^{k_z} e^{\alpha_k r_k^2} (z - Z_l)^{l_z+1} e^{\alpha_l r_l^2} 
                        + Z_l \int dxdydz (z - Z_k)^{k_z} e^{\alpha_k r_k^2} (z - Z_l)^{l_z} e^{\alpha_l r_l^2} 
                      = ao_part_1 + ao_part_2

with :math:`r_l = \sqrt{(x-X_l)^2 + (y-Y_l)^2 + (z - Z_l)^2)}`
'''

in_fid = 'h2o.molden'
# Read the input file
qc = main_read(in_fid,itype='molden',all_mo=False)

# Compute atomic orbital overlap matrix 
ao_overlap_matrix = get_ao_overlap(qc.geo_spec,qc.geo_spec,qc.ao_spec)

# Compute the overlap of the molecular orbitals and weight it with the occupation number
electron_number = 0.
for i_mo in qc.mo_spec:

  electron_number += i_mo['occ_num'] * get_mo_overlap(i_mo['coeffs'],