#!/usr/bin/env python3
"""
convert hdf5 output (e.g. from PySCF) to ezfio

Usage:
    qp_convert_h5_to_ezfio [--noqmc] [--rmg] [-o EZFIO_DIR] FILE

Options:
    -o --output=EZFIO_DIR    Produced directory
                             by default is FILE.ezfio
    --noqmc                  don't include basis, cell, etc. for QMCPACK
    --rmg                      h5 contains cholesky decomposition informatin, these h5 result from RMG and the pyscf AFQMC converter of QMCPACK.

"""
from ezfio import ezfio
import h5py
import sys
import numpy as np
import os
from docopt import docopt
import gzip
#fname = sys.argv[1]
#qph5name = sys.argv[2]
def idx2_tri(i,j):
    """
    for 0-indexed counting
    """
    p = max(i,j)
    q = min(i,j)
    return q + (p*(p+1))//2
def idx2_tri_1(i,j):
    """
    for 1-indexed counting
    """
    p = max(i,j)
    q = min(i,j)
    return q + (p*(p-1))//2

def idx4_cplx_1(i,j,k,l):
    """
    original function from qp2 (fortran counting)
    """
    p = idx2_tri_1(i,k)
    q = idx2_tri_1(j,l)
    i1 = idx2_tri_1(p,q)
    return (i1,p,q)

def ao_idx_map_sign(i,j,k,l):
    """
    qp2 indexing
    """
    idx,ik,jl = idx4_cplx_1(i,j,k,l)
    ij = idx2_tri_1(i,j)
    kl = idx2_tri_1(k,l)
    idx = 2*idx - 1

    if (ij==kl):
        sign = 0.0
        use_map1 = False
    else:
        if ik==jl:
            if i<k:
                sign = 1.0
                use_map1 = True
            else:
                sign = -1.0
                use_map1 = True
        elif i==k:
            if j<l:
                sign = 1.0
                use_map1 = True
            else:
                sign = -1.0
                use_map1 = True
        elif j==l:
            if i<k:
                sign = 1.0
                use_map1 = True
            else:
                sign = -1.0
                use_map1 = True
        elif ((i<k) == (j<l)):
            if i<k:
                sign = 1.0
                use_map1 = True
            else:
                sign = -1.0
                use_map1 = True
        else:
            if ((j<l) == (ik<jl)):
                sign = 1.0
                use_map1 = False
            else:
                sign = -1.0
                use_map1 = False
    return (idx, use_map1, sign)
def get_ao_int_cplx(i,j,k,l,map1,map2):
    idx,use_m1,sgn = ao_idx_map_sign(i,j,k,l)
    if use_m1:
        tmp_re = map1[idx]
        tmp_im = map1[idx+1]
        tmp_im *= sgn
    else:
        tmp_re = map2[idx]
        if sgn != 0.0:
            tmp_im = map2[idx+1]
            tmp_im *= sgn
        else:
            tmp_im = 0.0
    return tmp_re + 1j*tmp_im
def kconserv_p_from_qkk2_mk(qkk2,mk):
    nk, nk2 = qkk2.shape
    assert(nk == nk2)
    kcon_p = np.zeros((nk,nk,nk),dtype=int)
    for i in range(nk):
        for j in range(nk):
            for k in range(nk):
                kcon_p[i,j,k] = qkk2[mk[j],qkk2[k,i]]
                assert(qkk2[mk[j],qkk2[k,i]] == qkk2[qkk2[j,k],i])
    return kcon_p


def get_full_path(file_path):
    file_path = os.path.expanduser(file_path)
    file_path = os.path.expandvars(file_path)
#    file_path = os.path.abspath(file_path)
    return file_path

def make_reim_identity_kblocks(nk,nm,na=None):
    if na is None:
        na = nm
    single_block = np.eye(nm, na, dtype=np.complex128).view(dtype=np.float64).reshape((nm, na, 2))
    kblocks = np.tile(single_block,[nk, 1, 1, 1])
    return kblocks

def make_reim_identity_block_diag(nk,nm,na=None):
    from scipy.linalg import block_diag
    kblocks = make_reim_identity_kblocks(nk,nm,na).view(dtype=np.complex128).squeeze()
    kblockdiag = block_diag(*kblocks).view(dtype=np.float64).reshape((nk*nm,nk*na,2))
    print(f'kblockdiag.shape = {kblockdiag.shape}')
    return kblockdiag

def flatten(l):
    res = []
    for i in l:
        if hasattr(i, "__iter__") and not isinstance(i, str):
            res.extend(flatten(i))
        else:
            res.append(i)
    return res

def save_array_do(ezfioname,subdir,data,chunksize=16384):
    dims = list(reversed(data.shape))
    rank = len(dims)
    flatdata = data.reshape(-1)
    dim_max = 1
    for i in dims:
        dim_max *= i
    with gzip.open(os.path.join(ezfioname,subdir)+'.gz','wb') as f:
        f.write(f'{rank:3d}\n'.encode())
        for d in dims:
            f.write(f'{d:20d} '.encode())
        f.write("\n".encode())
        fmtstring = chunksize*'{:24.15E}\n'
        for i in range(dim_max//chunksize):
            #f.write((chunksize*'{:24.15E}\n').format(*flatdata[i*chunksize:(i+1)*chunksize]).encode())
            #f.write(fmtstring.format(*flatdata[i*chunksize:(i+1)*chunksize]).encode())
            f.write((''.join("%24.15E\n" % xi for xi in flatdata[i*chunksize:(i+1)*chunksize])).encode())
            print(f'{100.*i/(dim_max//chunksize):7.3f}% complete')
        rem = dim_max%chunksize
        if rem:
            f.write((rem*'{:24.15E}\n').format(*flatdata[-rem:]).encode())
    return


def convert_mol(filename,qph5path):
    ezfio.set_file(filename)
    ezfio.set_nuclei_is_complex(False)

    with h5py.File(qph5path,'r') as qph5:
        nucl_num = qph5['nuclei'].attrs['nucl_num']
        ao_num = qph5['ao_basis'].attrs['ao_num']
        mo_num = qph5['mo_basis'].attrs['mo_num']
        elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num']
        elec_beta_num = qph5['electrons'].attrs['elec_beta_num']

    ezfio.set_nuclei_nucl_num(nucl_num)

    ezfio.set_ao_basis_ao_num(ao_num)
    ezfio.set_mo_basis_mo_num(mo_num)
    ezfio.electrons_elec_alpha_num = elec_alpha_num
    ezfio.electrons_elec_beta_num  = elec_beta_num



    ##ao_num = mo_num
    ##Important !
    #import math
    #nelec_per_kpt = num_elec // n_kpts
    #nelec_alpha_per_kpt = int(math.ceil(nelec_per_kpt / 2.))
    #nelec_beta_per_kpt = int(math.floor(nelec_per_kpt / 2.))
    #
    #ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts)
    #ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts)

    #ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
    #ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))

    #ezfio.set_utils_num_kpts(n_kpts)
    #ezfio.set_integrals_bielec_df_num(n_aux)

    #(old)Important
    #ezfio.set_nuclei_nucl_num(nucl_num)
    #ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
    #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
    #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )


    with h5py.File(qph5path,'r') as qph5:
        nucl_charge=qph5['nuclei/nucl_charge'][()].tolist()
        nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist()
        nucl_label=qph5['nuclei/nucl_label'][()].tolist()
        nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion']

    ezfio.set_nuclei_nucl_charge(nucl_charge)
    ezfio.set_nuclei_nucl_coord(nucl_coord)
    if isinstance(nucl_label[0],bytes):
        nucl_label = list(map(lambda x:x.decode(),nucl_label))
    ezfio.set_nuclei_nucl_label(nucl_label)

    ezfio.set_nuclei_io_nuclear_repulsion('Read')
    ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion)


    ##########################################
    #                                        #
    #                Basis                   #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        do_pseudo = qph5['pseudo'].attrs['do_pseudo']
        ezfio.set_pseudo_do_pseudo(do_pseudo)
        if (do_pseudo):
            ezfio.set_pseudo_pseudo_lmax(qph5['pseudo'].attrs['pseudo_lmax'])
            ezfio.set_pseudo_pseudo_klocmax(qph5['pseudo'].attrs['pseudo_klocmax'])
            ezfio.set_pseudo_pseudo_kmax(qph5['pseudo'].attrs['pseudo_kmax'])
            ezfio.set_pseudo_nucl_charge_remove(qph5['pseudo/nucl_charge_remove'][()].tolist())
            ezfio.set_pseudo_pseudo_n_k(qph5['pseudo/pseudo_n_k'][()].tolist())
            ezfio.set_pseudo_pseudo_n_kl(qph5['pseudo/pseudo_n_kl'][()].tolist())
            ezfio.set_pseudo_pseudo_v_k(qph5['pseudo/pseudo_v_k'][()].tolist())
            ezfio.set_pseudo_pseudo_v_kl(qph5['pseudo/pseudo_v_kl'][()].tolist())
            ezfio.set_pseudo_pseudo_dz_k(qph5['pseudo/pseudo_dz_k'][()].tolist())
            ezfio.set_pseudo_pseudo_dz_kl(qph5['pseudo/pseudo_dz_kl'][()].tolist())

    ##########################################
    #                                        #
    #                Basis                   #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        #coeftmp = qph5['ao_basis/ao_coef'][()]
        #expotmp = qph5['ao_basis/ao_expo'][()]
        ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis'])
        ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist())
        ezfio.set_ao_basis_ao_prim_num(qph5['ao_basis/ao_prim_num'][()].tolist())
        ezfio.set_ao_basis_ao_power(qph5['ao_basis/ao_power'][()].tolist())
        ezfio.set_ao_basis_ao_coef(qph5['ao_basis/ao_coef'][()].tolist())
        ezfio.set_ao_basis_ao_expo(qph5['ao_basis/ao_expo'][()].tolist())


    ##########################################
    #                                        #
    #               MO Coef                  #
    #                                        #
    ##########################################


    with h5py.File(qph5path,'r') as qph5:
        mo_coef = qph5['mo_basis/mo_coef'][()].tolist()
    ezfio.set_mo_basis_mo_coef(mo_coef)
    #maybe fix qp so we don't need this?
    #ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num])

    return

def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
    import json
    from scipy.linalg import block_diag
    dump_fci, dump_cd, dump_fci2 = (False, False, False)

    ezfio.set_file(filename)
    ezfio.set_nuclei_is_complex(True)
    # Dummy atom since AFQMC h5 has no atom information
    #nucl_num = 1
    #ezfio.set_nuclei_nucl_num(nucl_num)
    #ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
    #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
    #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )
    with h5py.File(qph5path,'r') as qph5:
        kpt_num = qph5['Hamiltonian/KPoints'][()].shape[0]
        ham_dims = qph5['Hamiltonian/dims'][()]
        NMOPerKP = qph5['Hamiltonian/NMOPerKP'][()]
        _, _, kpt_num, orb_num, elec_alpha_num_tot, elec_beta_num_tot, _, nchol_maybe = ham_dims

        #for now, all kpts must have same number of MOs
        for nmoi in NMOPerKP:
            if nmoi != NMOPerKP[0]:
                print("ERROR: all KPs must have same number of MOs")
                raise ValueError

        #TODO: fix na, nb in rmg
        assert(elec_alpha_num_tot % kpt_num == 0)
        assert(elec_beta_num_tot % kpt_num == 0)
        elec_alpha_num_per_kpt = elec_alpha_num_tot // kpt_num
        elec_beta_num_per_kpt = elec_beta_num_tot // kpt_num
        #elec_alpha_num_per_kpt =  qph5['Hamiltonian']['dims'][4]
        #elec_beta_num_per_kpt =  qph5['Hamiltonian']['dims'][5]
        #orb_num =  qph5['Hamiltonian']['dims'][3]
        #try:
        #  is_ao =  json.loads(qph5['metadata'][()].decode("utf-8").replace("'",'"'))['ortho_ao']
        #  if is_ao:
        #    ao_num = orb_num
        #  elif is_ao ==False:
        #    mo_num = orb_num
        #  else:
        #      raise ValueError('Problem with ortho_ao key in metadata')
        #except:
        #    raise UnicodeDecodeError('metadata not correctly parsed from HDF5 file')


    ezfio.set_nuclei_kpt_num(kpt_num)
    kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2
    ezfio.set_nuclei_kpt_pair_num(kpt_pair_num)
    # don't multiply nuclei by kpt_num
    # work in k-space, not in equivalent supercell
    #nucl_num_per_kpt = nucl_num
    #ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)
    # these are totals (kpt_num * num_per_kpt)
    # need to change if we want to truncate orbital space within pyscf
    #if is_ao:
    #    ao_num = orb_num*kpt_num
    #TODO: fix this?
    ao_num_tot = orb_num
    ao_num_per_kpt = NMOPerKP[0]
    mo_num_tot = orb_num
    mo_num_per_kpt = NMOPerKP[0]
    #mo_num_per_kpt = ao_num_per_kpt
    ezfio.set_ao_basis_ao_num(ao_num_per_kpt * kpt_num)
    ezfio.set_mo_basis_mo_num(mo_num_per_kpt * kpt_num)
    ezfio.set_ao_basis_ao_num_per_kpt(ao_num_per_kpt)
    ezfio.set_mo_basis_mo_num_per_kpt(mo_num_per_kpt)
    ezfio.electrons_elec_alpha_num = elec_alpha_num_per_kpt * kpt_num
    ezfio.electrons_elec_beta_num  = elec_beta_num_per_kpt * kpt_num
    ##########################################
    #                                        #
    #                Basis                   #
    #                                        #
    ##########################################
    #TODO
    nucl_num = 1
    ezfio.set_nuclei_nucl_num(nucl_num)
    ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
    ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
    ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )
    nucl_num_per_kpt = nucl_num
    ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)
    ezfio.set_nuclei_io_kpt_symm('Read')
    ezfio.set_ao_basis_ao_basis("dummy basis")

    #nucleus on which each AO is centered
    ao_nucl = [1 for i in range(ao_num_per_kpt)]*kpt_num
    ezfio.set_ao_basis_ao_nucl(ao_nucl)


    #Just need one (can clean this up later)
    ao_prim_num_max  = 5

    d = [ [0] *ao_prim_num_max]*ao_num_tot
    ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num_tot)
    ezfio.set_ao_basis_ao_power(d)
    ezfio.set_ao_basis_ao_coef(d)
    ezfio.set_ao_basis_ao_expo(d)
    ##########################################
    #                                        #
    #               MOCoeff                  #
    #                                        #
    ##########################################
    #TODO
    #coef_per_kpt = np.eye(mo_num_per_kpt, ao_num_per_kpt, dtype=np.complex128).view(dtype=np.float64).reshape((mo_num_per_kpt, ao_num_per_kpt, 2))

    #mo_coef_kpts = np.tile(coef_per_kpt,[kpt_num, 1, 1, 1])
    #qph5.create_dataset('mo_basis/mo_coef_kpts',data=make_reim_identity_kblocks(kpt_num, mo_num_per_kpt, ao_num_per_kpt))
    ezfio.set_mo_basis_mo_coef_kpts(make_reim_identity_kblocks(kpt_num, mo_num_per_kpt, ao_num_per_kpt))
    ezfio.set_mo_basis_mo_coef_complex(make_reim_identity_block_diag(kpt_num, mo_num_per_kpt, ao_num_per_kpt))

    ##########################################
    #                                        #
    #            Integrals Mono              #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        ovlp_ao_reim = make_reim_identity_kblocks(kpt_num,ao_num_per_kpt,ao_num_per_kpt)
        ne_ao_reim = np.zeros((kpt_num,ao_num_per_kpt,ao_num_per_kpt,2),dtype=np.float64)
        kin_ao_reim = np.zeros((kpt_num,ao_num_per_kpt,ao_num_per_kpt,2),dtype=np.float64)
        try:
            #TODO: ensure pseudo is included in nuc-elec term?
            for i in range(kpt_num):
                ne_ao_reim[i] = qph5[f'Hamiltonian/H1_nuc_kp{i}'][()]
                kin_ao_reim[i] = qph5[f'Hamiltonian/H1_kinetic_kp{i}'][()]
        except KeyError:
            # QE/pyscf doesn't have separate kinetic, nuc-elec, pseudo 1e ints, so just combine in nuc-elec and set rest to zero
            for i in range(kpt_num):
                ne_ao_reim[i] = qph5[f'Hamiltonian/H1_kp{i}'][()]

        ezfio.set_ao_one_e_ints_ao_integrals_kinetic_kpts(kin_ao_reim)
        ezfio.set_ao_one_e_ints_ao_integrals_overlap_kpts(ovlp_ao_reim)
        ezfio.set_ao_one_e_ints_ao_integrals_n_e_kpts(ne_ao_reim)

        ezfio.set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
        ezfio.set_ao_one_e_ints_io_ao_integrals_overlap('Read')
        ezfio.set_ao_one_e_ints_io_ao_integrals_n_e('Read')


    """
    with h5py.File(qph5path,'r') as qph5:
        if is_ao:
            kin_ao_reim=
            ovlp_ao_reim=
            ne_ao_reim=

            ezfio.set_ao_one_e_ints_ao_integrals_kinetic_kpts(kin_ao_reim)
            ezfio.set_ao_one_e_ints_ao_integrals_overlap_kpts(ovlp_ao_reim)
            ezfio.set_ao_one_e_ints_ao_integrals_n_e_kpts(ne_ao_reim)

            ezfio.set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
            ezfio.set_ao_one_e_ints_io_ao_integrals_overlap('Read')
            ezfio.set_ao_one_e_ints_io_ao_integrals_n_e('Read')
        else:
            kin_mo_reim=
            ovlp_mo_reim=
            ne_mo_reim=

            ezfio.set_mo_one_e_ints_mo_integrals_kinetic_kpts(kin_mo_reim)
            ezfio.set_mo_one_e_ints_mo_integrals_overlap_kpts(ovlp_mo_reim)
            ezfio.set_mo_one_e_ints_mo_integrals_n_e_kpts(ne_mo_reim)

            ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read')
            ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read')
            ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read')
      """
    ##########################################
    #                                        #
    #               k-points                 #
    #                                        #
    ##########################################
    #TODO
    with h5py.File(qph5path,'r') as qph5:
        #kconserv = qph5['nuclei/kconserv'][()].tolist()
        minusk = qph5['Hamiltonian']['MinusK'][:]+1
        QKTok2 = qph5['Hamiltonian']['QKTok2'][:]+1
        #TODO: change this after rmg is fixed
        #minusk = QKTok2[:,0]
        kconserv = kconserv_p_from_qkk2_mk(QKTok2-1,minusk-1)+1
        unique_kpt_num = len(qph5['Hamiltonian']['KPFactorized'])
        unique_k_idx = []
        for i in qph5['Hamiltonian']['KPFactorized'].keys():
            unique_k_idx.append(int(i[1:])+1)
        unique_k_idx.sort()
    kpt_sparse_map = np.zeros(kpt_num,dtype=int)
    isparse=0
    #TODO: make robust: this assumes that for each pair, the one with data has a lower index
    for i in range(kpt_num):
        if i+1 in unique_k_idx:
            kpt_sparse_map[i] = isparse+1
            isparse += 1
        else:
            kpt_sparse_map[i] = -kpt_sparse_map[minusk[i]-1]
    ezfio.set_nuclei_kconserv(kconserv.transpose(2,1,0))
    ezfio.set_nuclei_io_kconserv('Read')
    ezfio.set_nuclei_minusk(minusk)
    ezfio.set_nuclei_qktok2(np.transpose(QKTok2)) #transposed for correct col-major ordering
    ezfio.set_nuclei_kpt_sparse_map(kpt_sparse_map)
    ezfio.set_nuclei_unique_kpt_num(unique_kpt_num)
    # kpt_sparse_map
    # unique_kpt_num
    # io_kpt_symm
    ##########################################
    #                                        #
    #             Integrals Bi               #
    #                                        #
    ##########################################

    # should this be in ao_basis? ao_two_e_ints?
    with h5py.File(qph5path,'r') as qph5:
        nchol_per_kpt_all =  qph5['Hamiltonian']['NCholPerKP'][:]
        print(f'nchol_per_kpt_full = {nchol_per_kpt_all}')
        #nchol_per_kpt = nchol_per_kpt_all[nchol_per_kpt_all != 0]
        nchol_per_kpt = nchol_per_kpt_all[np.array(unique_k_idx,dtype=int)-1]
        print(f'nchol_per_kpt = {nchol_per_kpt}')
        print(f'unique_k_idx = {unique_k_idx}')
        #for i in range(kpt_num):
        #    if i+1 in unique_k_idx:
        #        print('* ',i,nchol_per_kpt_all[i])
        #    else:
        #        print('  ',i,nchol_per_kpt_all[i])

        nchol_per_kpt_max = max(nchol_per_kpt)
        ezfio.set_ao_two_e_ints_chol_num(nchol_per_kpt)
        ezfio.set_ao_two_e_ints_chol_num_max(nchol_per_kpt_max)
        if is_ao:
            #ao_num_per_kpt = ao_num//kpt_num
            ezfio.set_ao_two_e_ints_io_chol_ao_integrals('Read')
            #ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt)))
            L_list = []
            L_all = np.zeros((unique_kpt_num, kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max,2),dtype=np.float64)
            print(f'kpt_sparse_map = {kpt_sparse_map}')
            print(f'unique_k_idx-1 = {np.array(unique_k_idx)-1}')
            for i in range(unique_kpt_num):
                ki = unique_k_idx[i]-1
                #print(i, ki)
                L_i = qph5[f'Hamiltonian/KPFactorized/L{ki}'][()].reshape((kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[kpt_sparse_map[ki]-1], 2))
                #L.reshape(kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[i], 2)
                #L = np.einsum("ijklm->ilkjm", A, B)
                L_all[i,:,:,:,:nchol_per_kpt[kpt_sparse_map[ki]-1],:] = L_i

                #(6, 5184, 2)
                """
                for cmplx in range(2):
                    for ao_idx_i in range(ao_num_per_kpt):
                        for ao_idx_j in range(ao_num_per_kpt):
                            for chol_idx in range(nchol_per_kpt[i]):
                                for kpt_idx in range(kpt_num):
                                    ao_chol_two_e_ints[cmplx][ao_idx_i][ao_idx_j][chol_idx][kpt_idx][i] = L[kpt_idx][ao_idx_i][ao_idx_j][chol_idx][cmplx]
                """
            #ao_chol_two_e_ints = np.vstack(L_list)
            #ao_chol_two_e_ints = ao_chol_two_e_ints.transpose()
            #TODO: check dims/reshape/transpose
            #ezfio.set_ao_two_e_ints_chol_ao_integrals_complex(L_all.transpose((5,2,3,4,1,0)))
            ezfio.set_ao_two_e_ints_chol_ao_integrals_complex(L_all.transpose((0,1,4,3,2,5)))

            def fortformat(x0):
                x = f'{abs(x0):25.14E}'.strip()
                xsign = '-' if (x0<0) else ''
                e = x.find('E')
                return xsign + f'0.{x[0]}{x[2:e]}{x[e:e+2]}{abs(int(x[e+1:])*1+1):02d}'
            if dump_cd:
                #for qi in range(unique_kpt_num):
                #    for ki in range(kpt_num):
                #        for i in range(ao_num_per_kpt):
                #            for j in range(ao_num_per_kpt):
                #                for ci in range(nchol_per_kpt_max):
                #                    vr = L_all[qi,ki,i,j,ci,0]
                #                    vi = L_all[qi,ki,i,j,ci,1]
                #                    print(f'{qi:6d} {ki:6d} {i:6d} {j:6d} {ci:6d} {vr:25.15E} {vi:25.15E}')
                Lnew = L_all.transpose((5,2,3,4,1,0))
                #for qi in range(unique_kpt_num):
                #    for ki in range(kpt_num):
                #        for ci in range(nchol_per_kpt_max):
                #            for j in range(ao_num_per_kpt):
                #                for i in range(ao_num_per_kpt):
                #                    vr = Lnew[0,i,j,ci,ki,qi]
                #                    vi = Lnew[1,i,j,ci,ki,qi]
                #                    print(f'{qi:6d} {ki:6d} {ci:6d} {j:6d} {i:6d} {vr:25.15E} {vi:25.15E}')
                for qi in range(unique_kpt_num):
                    for ki in range(kpt_num):
                        for ci in range(nchol_per_kpt_max):
                            for i in range(ao_num_per_kpt):
                                for j in range(ao_num_per_kpt):
                                    #vr = Lnew[0,i,j,ci,ki,qi]
                                    #vi = Lnew[1,i,j,ci,ki,qi]
                                    vr = L_all[qi,ki,i,j,ci,0]
                                    vi = L_all[qi,ki,i,j,ci,1]
                                    print(f'{qi+1:6d} {ki+1:6d} {ci+1:6d} {i+1:6d} {j+1:6d} {fortformat(vr):>25s} {fortformat(vi):>25s}')

            if dump_fci:
                Wfull = np.zeros((ao_num_tot, ao_num_tot, ao_num_tot, ao_num_tot), dtype=np.complex128)
                for Qi in range(kpt_num):
                    Qloc = abs(kpt_sparse_map[Qi])-1
                    Qneg = (kpt_sparse_map[Qi] < 0)
                    LQ00 = L_all[Qloc]
                    #LQ0a = LQ00.view(dtype=np.complex128)
                    #print(f'LQ0a.shape {LQ0a.shape}')
                    #LQ0a1 = LQ0a.reshape((kpt_num,ao_num_per_kpt,ao_num_per_kpt,-1))
                    #print(f'LQ0a1.shape {LQ0a1.shape}')
                    LQ0 = LQ00[:,:,:,:,0] + 1j*LQ00[:,:,:,:,1]
                    #print(f'LQ0.shape {LQ0.shape}')
                    #print(f'abdiff {np.abs(LQ0a1 - LQ0).max()}')

                    

                    for kp in range(kpt_num):
                        kr = QKTok2[Qi,kp]-1
                        for ks in range(kpt_num):
                            kq = QKTok2[Qi,ks]-1
                            # 3
                            #if Qneg:
                            #    A = LQ0[kr].transpose((1,0,2)).conj()
                            #    B = LQ0[kq]
                            #    W = np.einsum('prn,sqn->pqrs',A,B)
                            #else:
                            #    A = LQ0[kp]
                            #    B = LQ0[ks].transpose((1,0,2)).conj()
                            #    W = np.einsum('rpn,qsn->pqrs',A,B)
                            # 4
                            #if Qneg:
                            #    A = LQ0[kr].transpose((1,0,2)).conj()
                            #    B = LQ0[kq].transpose((1,0,2))
                            #    W = np.einsum('prn,sqn->pqrs',A,B)
                            #else:
                            #    A = LQ0[kp]
                            #    B = LQ0[ks].conj()
                            #    W = np.einsum('prn,sqn->pqrs',A,B)
                            # 5
                            #if Qneg:
                            #    A = LQ0[kr].transpose((1,0,2)).conj()
                            #    B = LQ0[kq].transpose((1,0,2))
                            #    W = np.einsum('prn,qsn->pqrs',A,B)
                            #else:
                            #    A = LQ0[kp]
                            #    B = LQ0[ks].conj()
                            #    W = np.einsum('prn,qsn->pqrs',A,B)
                            # 6
                            if Qneg:
                                A = LQ0[kr].transpose((1,0,2)).conj()
                                B = LQ0[kq].transpose((1,0,2))
                                W = np.einsum('prn,sqn->pqrs',A,B)
                            else:
                                A = LQ0[kp]
                                B = LQ0[ks].conj()
                                W = np.einsum('prn,sqn->pqrs',A,B)
                            p0 = kp*ao_num_per_kpt
                            r0 = kr*ao_num_per_kpt
                            q0 = kq*ao_num_per_kpt
                            s0 = ks*ao_num_per_kpt
                            for ip in range(ao_num_per_kpt):
                                for iq in range(ao_num_per_kpt):
                                    for ir in range(ao_num_per_kpt):
                                        for i_s in range(ao_num_per_kpt):
                                            v = W[ip,iq,ir,i_s]
                                            #print(f'{p0+ip:5d} {q0+iq:5d} {r0+ir:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}')
                                            print(f'{p0+ip:5d} {r0+ir:5d} {q0+iq:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}')
                            Wfull[p0:p0+ao_num_per_kpt,q0:q0+ao_num_per_kpt,r0:r0+ao_num_per_kpt,s0:s0+ao_num_per_kpt] = W.copy()
                H1 = np.zeros((ao_num_tot, ao_num_tot), dtype=np.complex128)
                for Qi in range(kpt_num):
                    hi0 = qph5[f'Hamiltonian/H1_kp{Qi}'][()]
                    hi = hi0[:,:,0] + 1j*hi0[:,:,1]
                    H1[Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt,Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt] = hi.copy()
                mo_occ = ([1,]* elec_alpha_num_per_kpt + [0,] * (ao_num_per_kpt - elec_alpha_num_per_kpt)) * kpt_num
                E1=0
                E2j=0
                E2k=0
                print("Jij Kij")
                for i in range(ao_num_tot):
                    for j in range(ao_num_tot):
                        print(f'{i:5d} {j:5d} {i:5d} {j:5d}  {Wfull[i,j,i,j].real:25.15E} {Wfull[i,j,i,j].imag:25.15E}')
                        print(f'{i:5d} {j:5d} {j:5d} {i:5d}  {Wfull[i,j,j,i].real:25.15E} {Wfull[i,j,j,i].imag:25.15E}')
                for imo, iocc in enumerate(mo_occ):
                    if iocc:
                        E1 += 2*H1[imo,imo]
                        for jmo, jocc in enumerate(mo_occ):
                            if jocc:
                                E2j += 2* Wfull[imo,jmo,imo,jmo]
                                E2k -= Wfull[imo,jmo,jmo,imo]
                print(f'E1  = {E1:25.15E}')
                print(f'E2j = {E2j:25.15E}')
                print(f'E2k = {E2k:25.15E}')
                print(f'E2  = {E2j+E2k:25.15E}')


            if dump_fci2:
                ao_map1 = {}
                ao_map2 = {}
                ao_map1_idx = {}
                ao_map2_idx = {}
                Wdict = {}
                for kQ in range(kpt_num):
                    for kl in range(kpt_num):
                        kj = QKTok2[kQ,kl]-1
                        if (kj>kl):
                            continue
                        kjkl2 = idx2_tri(kj,kl)
                        Qneg = (kpt_sparse_map[kQ] < 0)
                        Qloc = abs(kpt_sparse_map[kQ]) - 1
                        if not Qneg:
                            ints_jl0 = L_all[Qloc,kl,:,:,:,:]
                            ints_jl = ints_jl0[:,:,:,0]+1j*ints_jl0[:,:,:,1]
                        else:
                            ints_jl0 = L_all[Qloc,kj,:,:,:,:]
                            ints_jl1 = ints_jl0[:,:,:,0]+1j*ints_jl0[:,:,:,1]
                            ints_jl = ints_jl1.transpose((1,0,2)).conj()
                        for kk in range(kl+1):
                            ki = QKTok2[minusk[kk]-1,kQ]-1 #TODO: check
                            if ki != kconserv[kl,kk,kj]-1:
                                print(ki,kconserv[kl,kk,kj],kl,kk,kj)
                                assert( ki == kconserv[kl,kk,kj]-1) #TODO: check
                            if (ki > kl):
                                continue
                            kikk2 = idx2_tri(ki,kk)
                            if not Qneg:
                                ints_ik0 = L_all[Qloc,ki,:,:,:]
                                ints_ik = ints_ik0[:,:,:,0] + 1j*ints_ik0[:,:,:,1]
                            else:
                                ints_ik0 = L_all[Qloc,kk,:,:,:]
                                ints_ik1 = ints_ik0[:,:,:,0]+1j*ints_ik0[:,:,:,1]
                                ints_ik = ints_ik1.transpose((1,0,2)).conj()
                            ints_jl_flat = ints_jl.reshape((ao_num_per_kpt**2, -1))
                            ints_ik_flat = ints_ik.reshape((ao_num_per_kpt**2, -1))

                            ints_ikjl = np.einsum('an,bn->ab',ints_ik_flat,ints_jl_flat).reshape((ao_num_per_kpt,)*4)

                            for il in range(ao_num_per_kpt):
                                l = il + kl*ao_num_per_kpt
                                for ij in range(ao_num_per_kpt):
                                    j = ij + kj*ao_num_per_kpt
                                    if (j>l):
                                        break
                                    jl2 = idx2_tri(j,l)
                                    for ik in range(ao_num_per_kpt):
                                        k = ik + kk*ao_num_per_kpt
                                        if (k>l):
                                            break
                                        for ii in range(ao_num_per_kpt):
                                            i = ii + ki*ao_num_per_kpt
                                            if ((j==l) and (i>k)):
                                                break
                                            ik2 = idx2_tri(i,k)
                                            if (ik2 > jl2):
                                                break
                                            integral = ints_ikjl[ii,ik,ij,il]
                                            if abs(integral) < 1E-15:
                                                continue
                                            idx_tmp,use_map1,sign = ao_idx_map_sign(i+1,j+1,k+1,l+1)
                                            tmp_re = integral.real
                                            tmp_im = integral.imag
                                            Wdict[i,j,k,l] = tmp_re + 1j*tmp_im
                                            if use_map1:
                                                if idx_tmp in ao_map1:
                                                    print(idx_tmp,1)
                                                    raise
                                                ao_map1[idx_tmp] = tmp_re
                                                ao_map1_idx[idx_tmp] = (i,j,k,l,'re')
                                                if sign != 0.0:
                                                    ao_map1[idx_tmp+1] = tmp_im*sign
                                                    ao_map1_idx[idx_tmp+1] = (i,j,k,l,'im')
                                            else:
                                                if idx_tmp in ao_map2:
                                                    print(idx_tmp,2)
                                                    raise
                                                ao_map2[idx_tmp] = tmp_re
                                                ao_map2_idx[idx_tmp] = (i,j,k,l,'re')
                                                if sign != 0.0:
                                                    ao_map2[idx_tmp+1] = tmp_im*sign
                                                    ao_map2_idx[idx_tmp+1] = (i,j,k,l,'im')
#                for idx in ao_map1:
#                    i,j,k,l,ax = ao_map1_idx[idx]
#                    print(f'1,{idx},{i},{j},{k},{l},{ax},{ao_map1[idx]}')
#                for idx in ao_map2:
#                    i,j,k,l,ax = ao_map2_idx[idx]
#                    print(f'2,{idx},{i},{j},{k},{l},{ax},{ao_map2[idx]}')
                for idx in Wdict:
                    i,j,k,l = idx
                    v = Wdict[idx]
                    print(f'{i+1:6d} {j+1:6d} {k+1:6d} {l+1:6d} {fortformat(v.real):>25s} {fortformat(v.imag):>25s}')


                                            


                #for Qi in range(kpt_num):
                #    Qloc = abs(kpt_sparse_map[Qi])-1
                #    Qneg = (kpt_sparse_map[Qi] < 0)
                #    LQ00 = L_all[Qloc]
                #    #LQ0a = LQ00.view(dtype=np.complex128)
                #    #print(f'LQ0a.shape {LQ0a.shape}')
                #    #LQ0a1 = LQ0a.reshape((kpt_num,ao_num_per_kpt,ao_num_per_kpt,-1))
                #    #print(f'LQ0a1.shape {LQ0a1.shape}')
                #    LQ0 = LQ00[:,:,:,:,0] + 1j*LQ00[:,:,:,:,1]
                #    #print(f'LQ0.shape {LQ0.shape}')
                #    #print(f'abdiff {np.abs(LQ0a1 - LQ0).max()}')

                #    

                #    for kp in range(kpt_num):
                #        kr = QKTok2[Qi,kp]-1
                #        for ks in range(kpt_num):
                #            kq = QKTok2[Qi,ks]-1
                #            if Qneg:
                #                A = LQ0[kr].transpose((1,0,2)).conj()
                #                B = LQ0[kq]
                #                W = np.einsum('prn,sqn->pqrs',A,B)
                #            else:
                #                A = LQ0[kp]
                #                B = LQ0[ks].transpose((1,0,2)).conj()
                #                W = np.einsum('rpn,qsn->pqrs',A,B)
                #            p0 = kp*ao_num_per_kpt
                #            r0 = kr*ao_num_per_kpt
                #            q0 = kq*ao_num_per_kpt
                #            s0 = ks*ao_num_per_kpt
                #            for ip in range(ao_num_per_kpt):
                #                for iq in range(ao_num_per_kpt):
                #                    for ir in range(ao_num_per_kpt):
                #                        for i_s in range(ao_num_per_kpt):
                #                            v = W[ip,iq,ir,i_s]
                #                            print(f'{p0+ip:5d} {q0+iq:5d} {r0+ir:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}')
                #            Wfull[p0:p0+ao_num_per_kpt,q0:q0+ao_num_per_kpt,r0:r0+ao_num_per_kpt,s0:s0+ao_num_per_kpt] = W.copy()
                #H1 = np.zeros((ao_num_tot, ao_num_tot), dtype=np.complex128)
                #for Qi in range(kpt_num):
                #    hi0 = qph5[f'Hamiltonian/H1_kp{Qi}'][()]
                #    hi = hi0[:,:,0] + 1j*hi0[:,:,1]
                #    H1[Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt,Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt] = hi.copy()
                #mo_occ = ([1,]* elec_alpha_num_per_kpt + [0,] * (ao_num_per_kpt - elec_alpha_num_per_kpt)) * kpt_num
                #E1=0
                #E2j=0
                #E2k=0
                #for imo, iocc in enumerate(mo_occ):
                #    if iocc:
                #        E1 += 2*H1[imo,imo]
                #        for jmo, jocc in enumerate(mo_occ):
                #            if jocc:
                #                E2j += 2* Wfull[imo,jmo,imo,jmo]
                #                E2k -= Wfull[imo,jmo,jmo,imo]
                #print(f'E1  = {E1:25.15E}')
                #print(f'E2j = {E2j:25.15E}')
                #print(f'E2k = {E2k:25.15E}')
                #print(f'E2  = {E2j+E2k:25.15E}')


                            

            #(2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num#)
            """
            df_num = qph5['ao_two_e_ints'].attrs['df_num']
            ezfio.set_ao_two_e_ints_df_num(df_num)
            if 'df_ao_integrals' in qph5['ao_two_e_ints'].keys():
                dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist()
                ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim)
                ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read')
            """
        else:
            raise NotImplementedError
            """
            ezfio.set_io_chol_mo_integrals('Read')
            df_num = qph5['ao_two_e_ints'].attrs['df_num']
            ezfio.set_ao_two_e_ints_df_num(df_num)
            dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist()
            ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim)
            ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read')
            """
            #mo_num_per_kpt = ao_num//kpt_num
            ezfio.set_io_chol_mo_integrals('Read')
            #ao_chol_two_e_ints = np.zeros((2, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max, kpt_num, len(nchol_per_kpt)))
            L_list = []
            for i in len(nchol_per_kpt):
             L = qph5['Hamiltonian']['KPFactorized'][f'L{i}'][:]
             L.reshape(kpt_num, mo_num_per_kpt, mo_num_per_kpt, nchol_per_kpt[i], 2)
             L = np.einsum("ijklm->ilkjm", A, B)
             L_list.append(L)

             #(6, 5184, 2)
             """
             for cmplx in range(2):
                 for ao_idx_i in range(ao_num_per_kpt):
                     for ao_idx_j in range(ao_num_per_kpt):
                         for chol_idx in range(nchol_per_kpt[i]):
                             for kpt_idx in range(kpt_num):
                                 ao_chol_two_e_ints[cmplx][ao_idx_i][ao_idx_j][chol_idx][kpt_idx][i] = L[kpt_idx][ao_idx_i][ao_idx_j][chol_idx][cmplx]
             """
            mo_chol_two_e_ints = np.vstack(L_list)
            mo_chol_two_e_ints = mo_chol_two_e_ints.transpose()
            ezfio.set_chol_mo_integrals_complex(mo_chol_two_e_ints)
            return




def convert_kpts(filename,qph5path,qmcpack=True):
    ezfio.set_file(filename)
    ezfio.set_nuclei_is_complex(True)

    with h5py.File(qph5path,'r') as qph5:
        kpt_num = qph5['nuclei'].attrs['kpt_num']
        nucl_num = qph5['nuclei'].attrs['nucl_num']
        ao_num = qph5['ao_basis'].attrs['ao_num']
        mo_num = qph5['mo_basis'].attrs['mo_num']
        elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num']
        elec_beta_num = qph5['electrons'].attrs['elec_beta_num']

    ezfio.set_nuclei_kpt_num(kpt_num)
    kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2
    ezfio.set_nuclei_kpt_pair_num(kpt_pair_num)

    # don't multiply nuclei by kpt_num
    # work in k-space, not in equivalent supercell
    nucl_num_per_kpt = nucl_num
    ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)

    # these are totals (kpt_num * num_per_kpt)
    # need to change if we want to truncate orbital space within pyscf
    ezfio.set_ao_basis_ao_num(ao_num)
    ezfio.set_mo_basis_mo_num(mo_num)
    ezfio.set_ao_basis_ao_num_per_kpt(ao_num//kpt_num)
    ezfio.set_mo_basis_mo_num_per_kpt(mo_num//kpt_num)
    ezfio.electrons_elec_alpha_num = elec_alpha_num
    ezfio.electrons_elec_beta_num  = elec_beta_num



    ##ao_num = mo_num
    ##Important !
    #import math
    #nelec_per_kpt = num_elec // n_kpts
    #nelec_alpha_per_kpt = int(math.ceil(nelec_per_kpt / 2.))
    #nelec_beta_per_kpt = int(math.floor(nelec_per_kpt / 2.))
    #
    #ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts)
    #ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts)

    #ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
    #ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))

    #ezfio.set_utils_num_kpts(n_kpts)
    #ezfio.set_integrals_bielec_df_num(n_aux)

    #(old)Important
    #ezfio.set_nuclei_nucl_num(nucl_num)
    #ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
    #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
    #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )


    with h5py.File(qph5path,'r') as qph5:
        nucl_charge=qph5['nuclei/nucl_charge'][()].tolist()
        nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist()
        nucl_label=qph5['nuclei/nucl_label'][()].tolist()
        nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion']
        madconst = qph5['nuclei'].attrs.get('madelung_const',default=0)

    ezfio.set_nuclei_nucl_charge(nucl_charge)
    ezfio.set_nuclei_nucl_coord(nucl_coord)
    if isinstance(nucl_label[0],bytes):
        nucl_label = list(map(lambda x:x.decode(),nucl_label))
    ezfio.set_nuclei_nucl_label(nucl_label)

    ezfio.set_nuclei_io_nuclear_repulsion('Read')
    ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion)
    ezfio.set_nuclei_madelung_constant(madconst)


    ##########################################
    #                                        #
    #                Basis(Dummy)            #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis'])
        ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist())


    #Just need one (can clean this up later)
    ao_prim_num_max  = 5

    d = [ [0] *ao_prim_num_max]*ao_num
    ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num)
    ezfio.set_ao_basis_ao_power(d)
    ezfio.set_ao_basis_ao_coef(d)
    ezfio.set_ao_basis_ao_expo(d)


    ###########################################
    ##                                        #
    ##                Pseudo                  #
    ##                                        #
    ###########################################

    #with h5py.File(qph5path,'r') as qph5:
    #    do_pseudo = qph5['pseudo'].attrs['do_pseudo']
    #    ezfio.set_pseudo_do_pseudo(do_pseudo)
    #    if (do_pseudo):
    #        ezfio.set_pseudo_pseudo_lmax(qph5['pseudo'].attrs['pseudo_lmax'])
    #        ezfio.set_pseudo_pseudo_klocmax(qph5['pseudo'].attrs['pseudo_klocmax'])
    #        ezfio.set_pseudo_pseudo_kmax(qph5['pseudo'].attrs['pseudo_kmax'])
    #        ezfio.set_pseudo_nucl_charge_remove(qph5['pseudo/nucl_charge_remove'][()].tolist())
    #        ezfio.set_pseudo_pseudo_n_k(qph5['pseudo/pseudo_n_k'][()].tolist())
    #        ezfio.set_pseudo_pseudo_n_kl(qph5['pseudo/pseudo_n_kl'][()].tolist())
    #        ezfio.set_pseudo_pseudo_v_k(qph5['pseudo/pseudo_v_k'][()].tolist())
    #        ezfio.set_pseudo_pseudo_v_kl(qph5['pseudo/pseudo_v_kl'][()].tolist())
    #        ezfio.set_pseudo_pseudo_dz_k(qph5['pseudo/pseudo_dz_k'][()].tolist())
    #        ezfio.set_pseudo_pseudo_dz_kl(qph5['pseudo/pseudo_dz_kl'][()].tolist())

    ##########################################
    #                                        #
    #                Basis(Dummy)            #
    #                                        #
    ##########################################

    #with h5py.File(qph5path,'r') as qph5:
    #    coeftmp = qph5['ao_basis/ao_coef'][()]
    #    expotmp = qph5['ao_basis/ao_expo'][()]
    #    ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis'])
    #    ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist())
    #    ezfio.set_ao_basis_ao_prim_num(qph5['ao_basis/ao_prim_num'][()].tolist())
    #    ezfio.set_ao_basis_ao_power(qph5['ao_basis/ao_power'][()].tolist())
    #    ezfio.set_ao_basis_ao_coef(qph5['ao_basis/ao_coef'][()].tolist())
    #    ezfio.set_ao_basis_ao_expo(qph5['ao_basis/ao_expo'][()].tolist())


    ##########################################
    #                                        #
    #                Basis(QMC)              #
    #                                        #
    ##########################################

    if qmcpack:
        try:
            with h5py.File(qph5path,'r') as qph5:
                ezfio.set_qmcpack_qmc_nshell(qph5['qmcpack'].attrs['qmc_nshell'])
                ezfio.set_qmcpack_qmc_prim_num_max(qph5['qmcpack'].attrs['qmc_prim_num_max'])
                ezfio.set_qmcpack_qmc_nucl(qph5['qmcpack/qmc_nucl'][()].tolist())
                ezfio.set_qmcpack_qmc_prim_num(qph5['qmcpack/qmc_prim_num'][()].tolist())
                ezfio.set_qmcpack_qmc_lbas(qph5['qmcpack/qmc_lbas'][()].tolist())
                ezfio.set_qmcpack_qmc_coef(qph5['qmcpack/qmc_coef'][()].tolist())
                ezfio.set_qmcpack_qmc_expo(qph5['qmcpack/qmc_expo'][()].tolist())

                ezfio.set_qmcpack_qmc_pbc(qph5['qmcpack'].attrs['PBC'])
                ezfio.set_qmcpack_qmc_cart(qph5['qmcpack'].attrs['cart'])
                ezfio.set_qmcpack_qmc_pseudo(qph5['qmcpack'].attrs['Pseudo'])
                ezfio.set_qmcpack_supertwist(qph5['qmcpack/Super_Twist'][()].tolist())
                ezfio.set_qmcpack_latticevectors(qph5['qmcpack/LatticeVectors'][()].tolist())
                ezfio.set_qmcpack_qmc_phase(qph5['qmcpack/qmc_phase'][()].tolist())
                ezfio.set_qmcpack_qmc_mo_energy(qph5['qmcpack/eigenval'][()].tolist())

        except AttributeError as err:
            print("################################################")
            print("# ERROR: problem copying QMCPACK data to ezfio #")
            print("#        make sure qmcpack plugin is built     #")
            print("################################################")
            #print(f"AttributeError: {err}")
            print("to create ezfio without qmcpack data, use 'qp_convert_h5_to_ezfio --noqmc'")
            raise



    ##########################################
    #                                        #
    #               MO Coef                  #
    #                                        #
    ##########################################


    with h5py.File(qph5path,'r') as qph5:
        mo_coef_kpts = qph5['mo_basis/mo_coef_kpts'][()].tolist()
        mo_coef_cplx = qph5['mo_basis/mo_coef_complex'][()].tolist()
    ezfio.set_mo_basis_mo_coef_kpts(mo_coef_kpts)
    ezfio.set_mo_basis_mo_coef_complex(mo_coef_cplx)
    #maybe fix qp so we don't need this?
    #ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num])


    ##########################################
    #                                        #
    #            Integrals Mono              #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        if 'ao_one_e_ints' in qph5.keys():
            kin_ao_reim=qph5['ao_one_e_ints/ao_integrals_kinetic_kpts'][()].tolist()
            ovlp_ao_reim=qph5['ao_one_e_ints/ao_integrals_overlap_kpts'][()].tolist()
            ne_ao_reim=qph5['ao_one_e_ints/ao_integrals_n_e_kpts'][()].tolist()

            ezfio.set_ao_one_e_ints_ao_integrals_kinetic_kpts(kin_ao_reim)
            ezfio.set_ao_one_e_ints_ao_integrals_overlap_kpts(ovlp_ao_reim)
            ezfio.set_ao_one_e_ints_ao_integrals_n_e_kpts(ne_ao_reim)

            ezfio.set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
            ezfio.set_ao_one_e_ints_io_ao_integrals_overlap('Read')
            ezfio.set_ao_one_e_ints_io_ao_integrals_n_e('Read')


    with h5py.File(qph5path,'r') as qph5:
        if 'mo_one_e_ints' in qph5.keys():
            kin_mo_reim=qph5['mo_one_e_ints/mo_integrals_kinetic_kpts'][()].tolist()
            ovlp_mo_reim=qph5['mo_one_e_ints/mo_integrals_overlap_kpts'][()].tolist()
            ne_mo_reim=qph5['mo_one_e_ints/mo_integrals_n_e_kpts'][()].tolist()

            ezfio.set_mo_one_e_ints_mo_integrals_kinetic_kpts(kin_mo_reim)
            ezfio.set_mo_one_e_ints_mo_integrals_overlap_kpts(ovlp_mo_reim)
            #ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim)
            ezfio.set_mo_one_e_ints_mo_integrals_n_e_kpts(ne_mo_reim)

            ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read')
            ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read')
            #ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read')
            ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read')

    ##########################################
    #                                        #
    #               k-points                 #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        kconserv = qph5['nuclei/kconserv'][()].tolist()

    ezfio.set_nuclei_kconserv(kconserv)
    ezfio.set_nuclei_io_kconserv('Read')

    ##########################################
    #                                        #
    #             Integrals Bi               #
    #                                        #
    ##########################################

    # should this be in ao_basis? ao_two_e_ints?
    with h5py.File(qph5path,'r') as qph5:
        if 'ao_two_e_ints' in qph5.keys():
            df_num = qph5['ao_two_e_ints'].attrs['df_num']
            ezfio.set_ao_two_e_ints_df_num(df_num)
            if 'df_ao_integrals' in qph5['ao_two_e_ints'].keys():
    #            dfao_re0=qph5['ao_two_e_ints/df_ao_integrals_real'][()].transpose((3,2,1,0))
    #            dfao_im0=qph5['ao_two_e_ints/df_ao_integrals_imag'][()].transpose((3,2,1,0))
    #            dfao_cmplx0 = np.stack((dfao_re0,dfao_im0),axis=-1).tolist()
    #            ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_cmplx0)
                #dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist()
                dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()]
                save_array_do(filename,'ao_two_e_ints/df_ao_integrals_complex',dfao_reim)
                #ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim)
                #dfao_dims = list(reversed(dfao_reim.shape))
                #test_write_df_ao(,5,dfao_dims,dfao_reim.size,dfao_reim.ravel())

                ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read')

        if 'mo_two_e_ints' in qph5.keys():
            df_num = qph5['ao_two_e_ints'].attrs['df_num']
            ezfio.set_ao_two_e_ints_df_num(df_num)
    #        dfmo_re0=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0))
    #        dfmo_im0=qph5['mo_two_e_ints/df_mo_integrals_imag'][()].transpose((3,2,1,0))
    #        dfmo_cmplx0 = np.stack((dfmo_re0,dfmo_im0),axis=-1).tolist()
    #        ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_cmplx0)
            dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist()
            ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim)
            ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read')

    return

def convert_cplx(filename,qph5path):
    ezfio.set_file(filename)
    ezfio.set_nuclei_is_complex(True)

    with h5py.File(qph5path,'r') as qph5:
        kpt_num = qph5['nuclei'].attrs['kpt_num']
        nucl_num = qph5['nuclei'].attrs['nucl_num']
        ao_num = qph5['ao_basis'].attrs['ao_num']
        mo_num = qph5['mo_basis'].attrs['mo_num']
        elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num']
        elec_beta_num = qph5['electrons'].attrs['elec_beta_num']

    ezfio.set_nuclei_kpt_num(kpt_num)
    kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2
    ezfio.set_nuclei_kpt_pair_num(kpt_pair_num)

    # don't multiply nuclei by kpt_num
    # work in k-space, not in equivalent supercell
    nucl_num_per_kpt = nucl_num
    ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)

    # these are totals (kpt_num * num_per_kpt)
    # need to change if we want to truncate orbital space within pyscf
    ezfio.set_ao_basis_ao_num(ao_num)
    ezfio.set_mo_basis_mo_num(mo_num)
    ezfio.electrons_elec_alpha_num = elec_alpha_num
    ezfio.electrons_elec_beta_num  = elec_beta_num



    ##ao_num = mo_num
    ##Important !
    #import math
    #nelec_per_kpt = num_elec // n_kpts
    #nelec_alpha_per_kpt = int(math.ceil(nelec_per_kpt / 2.))
    #nelec_beta_per_kpt = int(math.floor(nelec_per_kpt / 2.))
    #
    #ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts)
    #ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts)

    #ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
    #ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))

    #ezfio.set_utils_num_kpts(n_kpts)
    #ezfio.set_integrals_bielec_df_num(n_aux)

    #(old)Important
    #ezfio.set_nuclei_nucl_num(nucl_num)
    #ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
    #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
    #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )


    with h5py.File(qph5path,'r') as qph5:
        nucl_charge=qph5['nuclei/nucl_charge'][()].tolist()
        nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist()
        nucl_label=qph5['nuclei/nucl_label'][()].tolist()
        nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion']

    ezfio.set_nuclei_nucl_charge(nucl_charge)
    ezfio.set_nuclei_nucl_coord(nucl_coord)
    if isinstance(nucl_label[0],bytes):
        nucl_label = list(map(lambda x:x.decode(),nucl_label))
    ezfio.set_nuclei_nucl_label(nucl_label)

    ezfio.set_nuclei_io_nuclear_repulsion('Read')
    ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion)


    ##########################################
    #                                        #
    #                Basis                   #
    #                                        #
    ##########################################

#    with h5py.File(qph5path,'r') as qph5:
#        ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis'])
#        ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist())
#
#
#    #Just need one (can clean this up later)
#    ao_prim_num_max  = 5
#
#    d = [ [0] *ao_prim_num_max]*ao_num
#    ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num)
#    ezfio.set_ao_basis_ao_power(d)
#    ezfio.set_ao_basis_ao_coef(d)
#    ezfio.set_ao_basis_ao_expo(d)

    ##########################################
    #                                        #
    #                Basis                   #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        do_pseudo = qph5['pseudo'].attrs['do_pseudo']
        ezfio.set_pseudo_do_pseudo(do_pseudo)
        if (do_pseudo):
            ezfio.set_pseudo_pseudo_lmax(qph5['pseudo'].attrs['pseudo_lmax'])
            ezfio.set_pseudo_pseudo_klocmax(qph5['pseudo'].attrs['pseudo_klocmax'])
            ezfio.set_pseudo_pseudo_kmax(qph5['pseudo'].attrs['pseudo_kmax'])
            ezfio.set_pseudo_nucl_charge_remove(qph5['pseudo/nucl_charge_remove'][()].tolist())
            ezfio.set_pseudo_pseudo_n_k(qph5['pseudo/pseudo_n_k'][()].tolist())
            ezfio.set_pseudo_pseudo_n_kl(qph5['pseudo/pseudo_n_kl'][()].tolist())
            ezfio.set_pseudo_pseudo_v_k(qph5['pseudo/pseudo_v_k'][()].tolist())
            ezfio.set_pseudo_pseudo_v_kl(qph5['pseudo/pseudo_v_kl'][()].tolist())
            ezfio.set_pseudo_pseudo_dz_k(qph5['pseudo/pseudo_dz_k'][()].tolist())
            ezfio.set_pseudo_pseudo_dz_kl(qph5['pseudo/pseudo_dz_kl'][()].tolist())

   ##########################################
   #                                        #
   #                Basis                   #
   #                                        #
   ##########################################

    with h5py.File(qph5path,'r') as qph5:
        coeftmp = qph5['ao_basis/ao_coef'][()]
        expotmp = qph5['ao_basis/ao_expo'][()]
        ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis'])
        ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist())
        ezfio.set_ao_basis_ao_prim_num(qph5['ao_basis/ao_prim_num'][()].tolist())
        ezfio.set_ao_basis_ao_power(qph5['ao_basis/ao_power'][()].tolist())
        ezfio.set_ao_basis_ao_coef(qph5['ao_basis/ao_coef'][()].tolist())
        ezfio.set_ao_basis_ao_expo(qph5['ao_basis/ao_expo'][()].tolist())




    ##########################################
    #                                        #
    #               MO Coef                  #
    #                                        #
    ##########################################


    with h5py.File(qph5path,'r') as qph5:
        mo_coef_reim = qph5['mo_basis/mo_coef_complex'][()].tolist()
    ezfio.set_mo_basis_mo_coef_complex(mo_coef_reim)
    #maybe fix qp so we don't need this?
    #ezfio.set_mo_basis_mo_coef([[i for i in range(mo_num)] * ao_num])


    ##########################################
    #                                        #
    #            Integrals Mono              #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        if 'ao_one_e_ints' in qph5.keys():
            kin_ao_reim=qph5['ao_one_e_ints/ao_integrals_kinetic'][()].tolist()
            ovlp_ao_reim=qph5['ao_one_e_ints/ao_integrals_overlap'][()].tolist()
            ne_ao_reim=qph5['ao_one_e_ints/ao_integrals_n_e'][()].tolist()

            ezfio.set_ao_one_e_ints_ao_integrals_kinetic_complex(kin_ao_reim)
            ezfio.set_ao_one_e_ints_ao_integrals_overlap_complex(ovlp_ao_reim)
            ezfio.set_ao_one_e_ints_ao_integrals_n_e_complex(ne_ao_reim)

            ezfio.set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
            ezfio.set_ao_one_e_ints_io_ao_integrals_overlap('Read')
            ezfio.set_ao_one_e_ints_io_ao_integrals_n_e('Read')


    with h5py.File(qph5path,'r') as qph5:
        if 'mo_one_e_ints' in qph5.keys():
            kin_mo_reim=qph5['mo_one_e_ints/mo_integrals_kinetic'][()].tolist()
            #ovlp_mo_reim=qph5['mo_one_e_ints/mo_integrals_overlap'][()].tolist()
            ne_mo_reim=qph5['mo_one_e_ints/mo_integrals_n_e'][()].tolist()

            ezfio.set_mo_one_e_ints_mo_integrals_kinetic_complex(kin_mo_reim)
            #ezfio.set_mo_one_e_ints_mo_integrals_overlap_complex(ovlp_mo_reim)
            #ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim)
            ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim)

            ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read')
            #ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read')
            #ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read')
            ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read')

    ##########################################
    #                                        #
    #               k-points                 #
    #                                        #
    ##########################################

    with h5py.File(qph5path,'r') as qph5:
        kconserv = qph5['nuclei/kconserv'][()].tolist()

    ezfio.set_nuclei_kconserv(kconserv)
    ezfio.set_nuclei_io_kconserv('Read')


    ##########################################
    #                                        #
    #             Integrals Bi               #
    #                                        #
    ##########################################

    # should this be in ao_basis? ao_two_e_ints?
    with h5py.File(qph5path,'r') as qph5:
        if 'ao_two_e_ints' in qph5.keys():
            df_num = qph5['ao_two_e_ints'].attrs['df_num']
            ezfio.set_ao_two_e_ints_df_num(df_num)
            if 'df_ao_integrals' in qph5['ao_two_e_ints'].keys():
    #            dfao_re0=qph5['ao_two_e_ints/df_ao_integrals_real'][()].transpose((3,2,1,0))
    #            dfao_im0=qph5['ao_two_e_ints/df_ao_integrals_imag'][()].transpose((3,2,1,0))
    #            dfao_cmplx0 = np.stack((dfao_re0,dfao_im0),axis=-1).tolist()
    #            ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_cmplx0)
                dfao_reim=qph5['ao_two_e_ints/df_ao_integrals'][()].tolist()
                ezfio.set_ao_two_e_ints_df_ao_integrals_complex(dfao_reim)
                ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read')

        if 'mo_two_e_ints' in qph5.keys():
            df_num = qph5['ao_two_e_ints'].attrs['df_num']
            ezfio.set_ao_two_e_ints_df_num(df_num)
    #        dfmo_re0=qph5['mo_two_e_ints/df_mo_integrals_real'][()].transpose((3,2,1,0))
    #        dfmo_im0=qph5['mo_two_e_ints/df_mo_integrals_imag'][()].transpose((3,2,1,0))
    #        dfmo_cmplx0 = np.stack((dfmo_re0,dfmo_im0),axis=-1).tolist()
    #        ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_cmplx0)
            dfmo_reim=qph5['mo_two_e_ints/df_mo_integrals'][()].tolist()
            ezfio.set_mo_two_e_ints_df_mo_integrals_complex(dfmo_reim)
            ezfio.set_mo_two_e_ints_io_df_mo_integrals('Read')

    return


if __name__ == '__main__':
    ARGUMENTS = docopt(__doc__)

    #for i in range(1,6):
    #    for j in range(1,6):
    #        for k in range(1,6):
    #            for l in range(1,6):
    #                idx,usem1,sgn = ao_idx_map_sign(i,j,k,l)
    #                print(f'{i:4d} {j:4d} {k:4d} {l:4d}      {str(usem1)[0]:s}  {idx:6d}     {sgn:5.1f}')
    FILE = get_full_path(ARGUMENTS['FILE'])
    qmcpack = True
    rmg = False
    if ARGUMENTS["--output"]:
        EZFIO_FILE = get_full_path(ARGUMENTS["--output"])
    else:
        EZFIO_FILE = "{0}.ezfio".format(FILE)

    if ARGUMENTS["--noqmc"]:
        qmcpack = False
    if ARGUMENTS["--rmg"]:
        rmg = True
    with h5py.File(FILE,'r') as qph5:
        try:
            do_kpts = ('kconserv' in qph5['nuclei'].keys())
        except:
            do_kpts = False
    if (do_kpts or rmg):
        print("converting HDF5 to EZFIO for periodic system")
        if rmg:
            print("Using RMG and AFQMC h5")
            convert_kpts_cd(EZFIO_FILE,FILE,qmcpack)
        else:
            convert_kpts(EZFIO_FILE,FILE,qmcpack)
    else:
        print("converting HDF5 to EZFIO for molecular system")
        convert_mol(EZFIO_FILE,FILE)

#    sys.stdout.flush()
#    if os.system("qp_run save_ortho_mos "+EZFIO_FILE) != 0:
#        print("""Warning: You need to run
#
#         qp run save_ortho_mos
#
#to be sure your MOs will be orthogonal, which is not the case when
#the MOs are read from output files (not enough precision in output).""")