diff --git a/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh b/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh index 9bce8816..94895f18 100755 --- a/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh +++ b/src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh @@ -17,7 +17,6 @@ echo 'Create EZFIO' qp_edit -c $ezfio &> /dev/null #cp $ezfio/{ao,mo}_basis/ao_md5 -qp_run import_kconserv $ezfio #qp_run import_ao_2e_complex $ezfio #qp_run dump_ao_2e_from_df $ezfio #Read the integral diff --git a/src/utils_complex/MolPyscfToQPkpts.py b/src/utils_complex/MolPyscfToQPkpts.py index f3119300..a9afbc27 100644 --- a/src/utils_complex/MolPyscfToQPkpts.py +++ b/src/utils_complex/MolPyscfToQPkpts.py @@ -681,6 +681,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, ########################################## with h5py.File(qph5path,'a') as qph5: + # k,mo,ao(,2) mo_coef_f = np.array(mo_k.transpose((0,2,1)),order='c') mo_coef_blocked=block_diag(*mo_k) mo_coef_blocked_f = block_diag(*mo_coef_f) @@ -688,8 +689,8 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, qph5.create_dataset('mo_basis/mo_coef_imag',data=mo_coef_blocked.imag) qph5.create_dataset('mo_basis/mo_coef_kpts_real',data=mo_k.real) qph5.create_dataset('mo_basis/mo_coef_kpts_imag',data=mo_k.imag) - qph5.create_dataset('mo_basis/mo_coef',data=mo_coef_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nao,2))) - qph5.create_dataset('mo_basis/mo_coef_kpts',data=mo_coef_f.view(dtype=np.float64).reshape((Nk,nmo,nao,2))) + qph5.create_dataset('mo_basis/mo_coef_complex',data=mo_coef_blocked_f.view(dtype=np.float64).reshape((Nk*nmo,Nk*nao,2))) + qph5.create_dataset('mo_basis/mo_coef_complex_kpts',data=mo_coef_f.view(dtype=np.float64).reshape((Nk,nmo,nao,2))) print_kpts_unblocked(mo_k,'C.qp',mo_coef_threshold) @@ -704,11 +705,20 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, ovlp_ao = get_ovlp_ao(mf) if print_ao_ints_mono: - kin_ao_blocked=block_diag(*kin_ao) - ovlp_ao_blocked=block_diag(*ovlp_ao) - ne_ao_blocked=block_diag(*ne_ao) with h5py.File(qph5path,'a') as qph5: + kin_ao_blocked=block_diag(*kin_ao) + ovlp_ao_blocked=block_diag(*ovlp_ao) + ne_ao_blocked=block_diag(*ne_ao) + + kin_ao_f = np.array(kin_ao.transpose((0,2,1)),order='c') + ovlp_ao_f = np.array(ovlp_ao.transpose((0,2,1)),order='c') + ne_ao_f = np.array(ne_ao.transpose((0,2,1)),order='c') + + kin_ao_blocked_f = block_diag(*kin_ao_f) + ovlp_ao_blocked_f = block_diag(*ovlp_ao_f) + ne_ao_blocked_f = block_diag(*ne_ao_f) + qph5.create_dataset('ao_one_e_ints/ao_integrals_kinetic_real',data=kin_ao_blocked.real) qph5.create_dataset('ao_one_e_ints/ao_integrals_kinetic_imag',data=kin_ao_blocked.imag) qph5.create_dataset('ao_one_e_ints/ao_integrals_overlap_real',data=ovlp_ao_blocked.real) @@ -716,6 +726,10 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, qph5.create_dataset('ao_one_e_ints/ao_integrals_n_e_real', data=ne_ao_blocked.real) qph5.create_dataset('ao_one_e_ints/ao_integrals_n_e_imag', data=ne_ao_blocked.imag) + qph5.create_dataset('ao_one_e_ints/ao_integrals_kinetic',data=kin_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2))) + qph5.create_dataset('ao_one_e_ints/ao_integrals_overlap',data=ovlp_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2))) + qph5.create_dataset('ao_one_e_ints/ao_integrals_n_e', data=ne_ao_blocked_f.view(dtype=np.float64).reshape((Nk*nao,Nk*nao,2))) + for fname,ints in zip(('S.qp','V.qp','T.qp'), (ovlp_ao, ne_ao, kin_ao)): print_kpts_unblocked_upper(ints,fname,thresh_mono) @@ -750,15 +764,8 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, kconserv = tools.get_kconserv(cell, kpts) with h5py.File(qph5path,'a') as qph5: - qph5.create_dataset('nuclei/kconserv',data=np.transpose(kconserv+1,(0,2,1))) - - kcon_test = np.zeros((Nk,Nk,Nk),dtype=int) - for a in range(Nk): - for b in range(Nk): - for c in range(Nk): - kcon_test[a,c,b] = kconserv[a,b,c]+1 - with h5py.File(qph5path,'a') as qph5: - qph5.create_dataset('nuclei/kconserv_test',data=kcon_test) + kcon_f_phys = np.array(kconserv.transpose((1,2,0)),order='c') + qph5.create_dataset('nuclei/kconserv',data=kcon_f_phys+1) print_kcon_chem_to_phys(kconserv,'K.qp') @@ -783,8 +790,8 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, j3ao_new = get_j3ao_new(mf.with_df._cderi,nao,Nk) with h5py.File(qph5path,'a') as qph5: - qph5.create_dataset('ao_two_e_ints/df_ao_integrals_real',data=j3arr.transpose((2,3,1,0)).real) - qph5.create_dataset('ao_two_e_ints/df_ao_integrals_imag',data=j3arr.transpose((2,3,1,0)).imag) + #qph5.create_dataset('ao_two_e_ints/df_ao_integrals_real',data=j3arr.transpose((2,3,1,0)).real) + #qph5.create_dataset('ao_two_e_ints/df_ao_integrals_imag',data=j3arr.transpose((2,3,1,0)).imag) qph5.create_dataset('ao_two_e_ints/df_ao_integrals',data=j3ao_new.view(dtype=np.float64).reshape((nkpt_pairs,naux,nao,nao,2))) if print_mo_ints_df: @@ -795,8 +802,8 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, print_df(j3mo,'D.mo.qp',bielec_int_threshold) with h5py.File(qph5path,'a') as qph5: - qph5.create_dataset('mo_two_e_ints/df_mo_integrals_real',data=j3mo.transpose((2,3,1,0)).real) - qph5.create_dataset('mo_two_e_ints/df_mo_integrals_imag',data=j3mo.transpose((2,3,1,0)).imag) + #qph5.create_dataset('mo_two_e_ints/df_mo_integrals_real',data=j3mo.transpose((2,3,1,0)).real) + #qph5.create_dataset('mo_two_e_ints/df_mo_integrals_imag',data=j3mo.transpose((2,3,1,0)).imag) qph5.create_dataset('mo_two_e_ints/df_mo_integrals',data=j3mo_new.view(dtype=np.float64).reshape((nkpt_pairs,naux,nmo,nmo,2))) if (print_ao_ints_bi): diff --git a/src/utils_complex/create_ezfio_complex_3idx.py b/src/utils_complex/create_ezfio_complex_3idx.py index a1a2cca9..34d2c801 100755 --- a/src/utils_complex/create_ezfio_complex_3idx.py +++ b/src/utils_complex/create_ezfio_complex_3idx.py @@ -5,33 +5,37 @@ import h5py import sys import numpy as np filename = sys.argv[1] -h5filename = sys.argv[2] -#num_elec, nucl_num, mo_num = map(int,sys.argv[2:5]) +qph5path = sys.argv[2] -#nuclear_repulsion = float(sys.argv[5]) -#ao_num = int(sys.argv[6]) -#n_kpts = int(sys.argv[7]) -#n_aux = int(sys.argv[8]) ezfio.set_file(filename) -qph5=h5py.File(h5filename,'r') +#qph5=h5py.File(qph5path,'r') -kpt_num = qph5['nuclei'].attrs['kpt_num'] +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.electrons_elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num'] -ezfio.electrons_elec_beta_num = qph5['electrons'].attrs['elec_beta_num'] -nucl_num = qph5['nuclei'].attrs['nucl_num'] -nucl_num_per_kpt = nucl_num // kpt_num -ao_num = qph5['ao_basis'].attrs['ao_num'] -mo_num = qph5['mo_basis'].attrs['mo_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 @@ -57,30 +61,30 @@ ezfio.set_mo_basis_mo_num(mo_num) #ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num ) #ezfio.set_nuclei_nucl_label( ['He'] * nucl_num ) -ezfio.set_nuclei_nucl_num(nucl_num_per_kpt) -nucl_charge=qph5['nuclei/nucl_charge'][()].tolist() +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) - -nucl_coord=qph5['nuclei/nucl_coord'][()].T.tolist() ezfio.set_nuclei_nucl_coord(nucl_coord) - -nucl_label=qph5['nuclei/nucl_label'][()].tolist() ezfio.set_nuclei_nucl_label(nucl_label) - ezfio.set_nuclei_io_nuclear_repulsion('Read') -nuclear_repulsion = qph5['nuclei'].attrs['nuclear_repulsion'] ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion) -# Ao num -#ao_num = mo_num -#ezfio.set_ao_basis_ao_basis("Dummy one. We read MO") -ezfio.set_ao_basis_ao_num(ao_num) -#ezfio.set_ao_basis_ao_nucl([1]*ao_num) #Maybe put a realy incorrect stuff -ezfio.set_ao_basis_ao_basis(qph5['ao_basis'].attrs['ao_basis']) -ezfio.set_ao_basis_ao_nucl(qph5['ao_basis/ao_nucl'][()].tolist()) +########################################## +# # +# 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) @@ -94,67 +98,82 @@ ezfio.set_ao_basis_ao_expo(d) -ezfio.set_mo_basis_mo_num(mo_num) -#c_mo = [[1 if i==j else 0 for i in range(mo_num)] for j in range(ao_num)] -#ezfio.set_mo_basis_mo_coef([ [0]*mo_num] * ao_num) -##ezfio.set_mo_basis_mo_coef_real(c_mo) - -mo_coef_re0 = qph5['mo_basis/mo_coef_real'][()].T -mo_coef_im0 = qph5['mo_basis/mo_coef_imag'][()].T -mo_coef_cmplx0 = np.stack((mo_coef_re0,mo_coef_im0),axis=-1).tolist() - -#ezfio.set_mo_basis_mo_coef_real(qph5['mo_basis/mo_coef_real'][()].tolist()) -#ezfio.set_mo_basis_mo_coef_imag(qph5['mo_basis/mo_coef_imag'][()].tolist()) -ezfio.set_mo_basis_mo_coef_complex(mo_coef_cmplx0) + +########################################## +# # +# 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]) -ezfio.set_nuclei_is_complex(True) -# fortran-ordered re,im parts -kin_ao_re0=qph5['ao_one_e_ints/ao_integrals_kinetic_real'][()].T -kin_ao_im0=qph5['ao_one_e_ints/ao_integrals_kinetic_imag'][()].T -#test where to stack? (axis=0 or -1?) -kin_ao_cmplx0=np.stack((kin_ao_re0,kin_ao_im0),axis=-1).tolist() +########################################## +# # +# Integrals Mono # +# # +########################################## + +with h5py.File(qph5path,'r') as qph5: + 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() -ovlp_ao_re0=qph5['ao_one_e_ints/ao_integrals_overlap_real'][()].T -ovlp_ao_im0=qph5['ao_one_e_ints/ao_integrals_overlap_imag'][()].T -#test where to stack? (axis=0 or -1?) -ovlp_ao_cmplx0=np.stack((ovlp_ao_re0,ovlp_ao_im0),axis=-1).tolist() - -ne_ao_re0=qph5['ao_one_e_ints/ao_integrals_n_e_real'][()].T -ne_ao_im0=qph5['ao_one_e_ints/ao_integrals_n_e_imag'][()].T -#test where to stack? (axis=0 or -1?) -ne_ao_cmplx0=np.stack((ne_ao_re0,ne_ao_im0),axis=-1).tolist() - -ezfio.set_ao_one_e_ints_ao_integrals_kinetic_complex(kin_ao_cmplx0) -ezfio.set_ao_one_e_ints_ao_integrals_overlap_complex(ovlp_ao_cmplx0) -ezfio.set_ao_one_e_ints_ao_integrals_n_e_complex(ne_ao_cmplx0) +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') + +########################################## +# # +# 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? -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_real' 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) - ezfio.set_ao_two_e_ints_io_df_ao_integrals('Read') +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) - ezfio.set_mo_two_e_ints_io_df_mo_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') #TODO: add check and only do this if ints exist diff --git a/src/utils_complex/dump_ao_1e_cplx.irp.f b/src/utils_complex/dump_ao_1e_cplx.irp.f index f49b2529..de5a48ee 100644 --- a/src/utils_complex/dump_ao_1e_cplx.irp.f +++ b/src/utils_complex/dump_ao_1e_cplx.irp.f @@ -13,6 +13,16 @@ subroutine run do i=1,ao_num write(*,'(200(E24.15))') ao_one_e_integrals_complex(i,:) enddo + write(*,'(A)') 'ao_kinetic_integrals_complex' + write(*,'(A)') '---------------' + do i=1,ao_num + write(*,'(200(E24.15))') ao_kinetic_integrals_complex(i,:) + enddo + write(*,'(A)') 'ao_ne_integrals_complex' + write(*,'(A)') '---------------' + do i=1,ao_num + write(*,'(200(E24.15))') ao_integrals_n_e_complex(i,:) + enddo write(*,'(A)') 'ao_overlap_complex' write(*,'(A)') '---------------' do i=1,ao_num