diff --git a/src/utils_complex/MolPyscfToQPkpts.py b/src/utils_complex/MolPyscfToQPkpts.py index a3a84235..71da37db 100644 --- a/src/utils_complex/MolPyscfToQPkpts.py +++ b/src/utils_complex/MolPyscfToQPkpts.py @@ -33,6 +33,9 @@ def idx40(i,j,k,l): def idx4(i,j,k,l): return idx2_tri((idx2_tri((i-1,k-1)),idx2_tri((j-1,l-1))))+1 +def stri4(i,j,k,l): + return (4*'{:5d}').format(i,j,k,l) + def stri4z(i,j,k,l,zr,zi): return (4*'{:5d}'+2*'{:25.16e}').format(i,j,k,l,zr,zi) @@ -509,6 +512,25 @@ def print_ao_bi(mf,kconserv=None,outfilename='W.ao.qp',bielec_int_threshold = 1E v.real,v.imag)+'\n') +def print_kcon_chem_to_phys(kcon,fname): + ''' + input: kconserv in chem notation kcon_c[a,b,c] = d + where (ab|cd) is allowed by symmetry + output: kconserv in phys notation kcon_p[i,j,k] = l + where is allowed by symmetry + (printed to file) + ''' + Nk,n2,n3 = kcon.shape + if (n2!=n3 or Nk!=n2): + raise Exception('print_kcon_chem_to_phys called with non-cubic array') + + with open(fname,'w') as outfile: + for a in range(Nk): + for b in range(Nk): + for c in range(Nk): + d = kcon[a,b,c] + outfile.write(stri4(a+1,c+1,b+1,d+1)+'\n') + def print_kpts_unblocked(ints_k,outfilename,thresh): ''' for ints_k of shape (Nk,n1,n2), @@ -556,26 +578,70 @@ def get_kin_ao(mf): def get_ovlp_ao(mf): nao = mf.cell.nao_nr() Nk = len(mf.kpts) - return np.reshape(mf.get_ovlp(cell=cell,kpts=kpts),(Nk,nao,nao)) + return np.reshape(mf.get_ovlp(cell=mf.cell,kpts=mf.kpts),(Nk,nao,nao)) def get_pot_ao(mf): nao = mf.cell.nao_nr() Nk = len(mf.kpts) if mf.cell.pseudo: - v_kpts_ao = np.reshape(mf.with_df.get_pp(kpts=kpts),(Nk,nao,nao)) + v_kpts_ao = np.reshape(mf.with_df.get_pp(kpts=mf.kpts),(Nk,nao,nao)) else: - v_kpts_ao = np.reshape(mf.with_df.get_nuc(kpts=kpts),(Nk,nao,nao)) + v_kpts_ao = np.reshape(mf.with_df.get_nuc(kpts=mf.kpts),(Nk,nao,nao)) - if len(cell._ecpbas) > 0: + if len(mf.cell._ecpbas) > 0: from pyscf.pbc.gto import ecp - v_kpts_ao += np.reshape(ecp.ecp_int(cell, kpts),(Nk,nao,nao)) + v_kpts_ao += np.reshape(ecp.ecp_int(mf.cell, mf.kpts),(Nk,nao,nao)) return v_kpts_ao def ao_to_mo_1e(ao_kpts,mo_coef): return np.einsum('kim,kij,kjn->kmn',mo_coef.conj(),ao_kpts_ao,mo_coef) +def get_j3ao(fname,nao,Nk): + import h5py + with h5py.File(fname,'r') as intfile: + j3c = intfile.get('j3c') + j3ckeys = list(j3c.keys()) + j3ckeys.sort(key=lambda strkey:int(strkey)) + + # in new(?) version of PySCF, there is an extra layer of groups before the datasets + # datasets used to be [/j3c/0, /j3c/1, /j3c/2, ...] + # datasets now are [/j3c/0/0, /j3c/1/0, /j3c/2/0, ...] + j3clist = [j3c.get(i+'/0') for i in j3ckeys] + #if j3clist==[None]*len(j3clist): + if not(any(j3clist)): + # if using older version, stop before last level + j3clist = [j3c.get(i) for i in j3ckeys] + + naosq = nao*nao + naotri = (nao*(nao+1))//2 + nkinvsq = 1./np.sqrt(Nk) + + # dimensions are (kikj,iaux,jao,kao), where kikj is compound index of kpts i and j + # output dimensions should be reversed (nao, nao, naux, nkptpairs) + return np.array([(i.value.reshape([-1,nao,nao]) if (i.shape[1] == naosq) else makesq3(i.value,nao)) * nkinvsq for i in j3clist]) + +def print_df(j3arr,fname,thresh): + with open(fname,'w') as outfile: + for k,kpt_pair in enumerate(j3arr): + for iaux,dfbasfunc in enumerate(kpt_pair): + for i,i0 in enumerate(dfbasfunc): + for j,v in enumerate(i0): + if (abs(v) > thresh): + outfile.write(stri4z(i+1,j+1,iaux+1,k+1,v.real,v.imag)+'\n') + return + +def df_pad_ref_test(j3arr,nao,naux,nkpt_pairs): + df_ao_tmp = np.zeros((nao,nao,naux,nkpt_pairs),dtype=np.complex128) + for k,kpt_pair in enumerate(j3arr): + for iaux,dfbasfunc in enumerate(kpt_pair): + for i,i0 in enumerate(dfbasfunc): + for j,v in enumerate(i0): + df_ao_tmp[i,j,iaux,k]=v + return df_ao_tmp + + def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, print_ao_ints_bi=False, @@ -591,12 +657,11 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, int_threshold = The integral will be not printed in they are bellow that ''' - from pyscf.pbc import ao2mo +# from pyscf.pbc import ao2mo from pyscf.pbc import tools - from pyscf.pbc.gto import ecp - from pyscf.data import nist import h5py - import scipy +# import scipy + from scipy.linalg import block_diag mo_coef_threshold = int_threshold ovlp_threshold = int_threshold @@ -614,98 +679,127 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, qph5.create_group('ao_basis') qph5.create_group('mo_basis') - qph5 = h5py.File(qph5path,'a') - natom = cell.natm - nelec = cell.nelectron - neleca,nelecb = cell.nelec - atom_xyz = mf.cell.atom_coords() - if not(mf.cell.unit.startswith(('B','b','au','AU'))): - atom_xyz /= nist.BOHR # always convert to au - - strtype=h5py.special_dtype(vlen=str) - atom_dset=qph5.create_dataset('nuclei/nucl_label',(natom,),dtype=strtype) - for i in range(natom): - atom_dset[i] = mf.cell.atom_pure_symbol(i) - qph5.create_dataset('nuclei/nucl_coord',data=atom_xyz) - qph5.create_dataset('nuclei/nucl_charge',data=mf.cell.atom_charges()) - - - print('n_atom per kpt', natom) - print('num_elec per kpt', nelec) - mo_coeff = mf.mo_coeff # Mo_coeff actif mo_k = np.array([c[:,cas_idx] for c in mo_coeff] if cas_idx is not None else mo_coeff) e_k = np.array([e[cas_idx] for e in mf.mo_energy] if cas_idx is not None else mf.mo_energy) Nk, nao, nmo = mo_k.shape + print("n Kpts", Nk) print("n active Mos per kpt", nmo) print("n AOs per kpt", nao) -# naux = mf.with_df.auxcell.nao -# print("n df fitting functions", naux) - - #in old version: param << nelec*Nk, nmo*Nk, natom*Nk - qph5['electrons'].attrs['elec_alpha_num']=neleca*Nk - qph5['electrons'].attrs['elec_beta_num']=nelecb*Nk - qph5['mo_basis'].attrs['mo_num']=Nk*nmo - qph5['ao_basis'].attrs['ao_num']=Nk*nao - #qph5['nuclei'].attrs['nucl_num']=Nk*natom - qph5['nuclei'].attrs['nucl_num']=natom - qph5['nuclei'].attrs['kpt_num']=Nk - qph5.create_group('ao_two_e_ints') -# qph5['ao_two_e_ints'].attrs['df_num']=naux + ########################################## + # # + # Nuclei # + # # + ########################################## - qph5['ao_basis'].attrs['ao_basis']=mf.cell.basis - ao_nucl=[mf.cell.bas_atom(i)+1 for i in range(nao)] - qph5.create_dataset('ao_basis/ao_nucl',data=Nk*ao_nucl) + natom = cell.natm + print('n_atom per kpt', natom) + atom_xyz = mf.cell.atom_coords() + if not(mf.cell.unit.startswith(('B','b','au','AU'))): + from pyscf.data.nist import BOHR + atom_xyz /= BOHR # always convert to au + + with h5py.File(qph5path,'a') as qph5: + qph5['nuclei'].attrs['kpt_num']=Nk + qph5['nuclei'].attrs['nucl_num']=natom + qph5.create_dataset('nuclei/nucl_coord',data=atom_xyz) + qph5.create_dataset('nuclei/nucl_charge',data=mf.cell.atom_charges()) + + strtype=h5py.special_dtype(vlen=str) + atom_dset=qph5.create_dataset('nuclei/nucl_label',(natom,),dtype=strtype) + for i in range(natom): + atom_dset[i] = mf.cell.atom_pure_symbol(i) + + ########################################## + # # + # Basis # + # # + ########################################## + + # nucleus on which each AO is centered + ao_nucl=[i[0] for i in mf.cell.ao_labels(fmt=False,base=1)] + + with h5py.File(qph5path,'a') as qph5: + qph5['mo_basis'].attrs['mo_num']=Nk*nmo + qph5['ao_basis'].attrs['ao_num']=Nk*nao + + qph5['ao_basis'].attrs['ao_basis']=mf.cell.basis + + qph5.create_dataset('ao_basis/ao_nucl',data=Nk*ao_nucl) + + ########################################## + # # + # Electrons # + # # + ########################################## + + nelec = cell.nelectron + neleca,nelecb = cell.nelec + + print('num_elec per kpt', nelec) + + with h5py.File(qph5path,'a') as qph5: + #in old version: param << nelec*Nk, nmo*Nk, natom*Nk + qph5['electrons'].attrs['elec_alpha_num']=neleca*Nk + qph5['electrons'].attrs['elec_beta_num']=nelecb*Nk + + ########################################## + # # + # Nuclear Repulsion # + # # + ########################################## - # _ - # |\ | _ | _ _. ._ |_) _ ._ | _ o _ ._ - # | \| |_| (_ | (/_ (_| | | \ (/_ |_) |_| | _> | (_) | | - # | - #Total energy shift due to Ewald probe charge = -1/2 * Nelec*madelung/cell.vol = shift = tools.pbc.madelung(cell, kpts)*cell.nelectron * -.5 e_nuc = (cell.energy_nuc() + shift)*Nk print('nucl_repul', e_nuc) - qph5['nuclei'].attrs['nuclear_repulsion']=e_nuc + + with h5py.File(qph5path,'a') as qph5: + qph5['nuclei'].attrs['nuclear_repulsion']=e_nuc - # __ __ _ - # |\/| | | | _ _ |_ _ - # | | |__| |__ (_) (/_ | _> - # - mo_coef_blocked=scipy.linalg.block_diag(*mo_k) - qph5.create_dataset('mo_basis/mo_coef_real',data=mo_coef_blocked.real) - 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) + ########################################## + # # + # MO Coef # + # # + ########################################## + + mo_coef_blocked=block_diag(*mo_k) + with h5py.File(qph5path,'a') as qph5: + qph5.create_dataset('mo_basis/mo_coef_real',data=mo_coef_blocked.real) + 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) print_kpts_unblocked(mo_k,'C.qp',mo_coef_threshold) - # ___ - # | ._ _|_ _ _ ._ _. | _ |\/| _ ._ _ - # _|_ | | |_ (/_ (_| | (_| | _> | | (_) | | (_) - # _| + ########################################## + # # + # Integrals Mono # + # # + ########################################## ne_ao = get_pot_ao(mf) kin_ao = get_kin_ao(mf) ovlp_ao = get_ovlp_ao(mf) if print_ao_ints_mono: - kin_ao_blocked=scipy.linalg.block_diag(*kin_ao) - ovlp_ao_blocked=scipy.linalg.block_diag(*ovlp_ao) - ne_ao_blocked=scipy.linalg.block_diag(*ne_ao) + kin_ao_blocked=block_diag(*kin_ao) + ovlp_ao_blocked=block_diag(*ovlp_ao) + ne_ao_blocked=block_diag(*ne_ao) - 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) - qph5.create_dataset('ao_one_e_ints/ao_integrals_overlap_imag',data=ovlp_ao_blocked.imag) - 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) + with h5py.File(qph5path,'a') as qph5: + 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) + qph5.create_dataset('ao_one_e_ints/ao_integrals_overlap_imag',data=ovlp_ao_blocked.imag) + 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) for fname,ints in zip(('S.qp','V.qp','T.qp'), (ovlp_ao, ne_ao, kin_ao)): @@ -716,107 +810,91 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, ovlp_mo = ao_to_mo_1e(ovlp_ao,mo_k) ne_mo = ao_to_mo_1e(ne_ao,mo_k) - kin_mo_blocked=scipy.linalg.block_diag(*kin_mo) - ovlp_mo_blocked=scipy.linalg.block_diag(*ovlp_mo) - ne_mo_blocked=scipy.linalg.block_diag(*ne_mo) + kin_mo_blocked=block_diag(*kin_mo) + ovlp_mo_blocked=block_diag(*ovlp_mo) + ne_mo_blocked=block_diag(*ne_mo) - qph5.create_dataset('mo_one_e_ints/mo_integrals_kinetic_real',data=kin_mo_blocked.real) - qph5.create_dataset('mo_one_e_ints/mo_integrals_kinetic_imag',data=kin_mo_blocked.imag) - qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap_real',data=ovlp_mo_blocked.real) - qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap_imag',data=ovlp_mo_blocked.imag) - qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_real', data=ne_mo_blocked.real) - qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_imag', data=ne_mo_blocked.imag) + with h5py.File(qph5path,'a') as qph5: + qph5.create_dataset('mo_one_e_ints/mo_integrals_kinetic_real',data=kin_mo_blocked.real) + qph5.create_dataset('mo_one_e_ints/mo_integrals_kinetic_imag',data=kin_mo_blocked.imag) + qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap_real',data=ovlp_mo_blocked.real) + qph5.create_dataset('mo_one_e_ints/mo_integrals_overlap_imag',data=ovlp_mo_blocked.imag) + qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_real', data=ne_mo_blocked.real) + qph5.create_dataset('mo_one_e_ints/mo_integrals_n_e_imag', data=ne_mo_blocked.imag) for fname,ints in zip(('S.mo.qp','V.mo.qp','T.mo.qp'), (ovlp_mo, ne_mo, kin_mo)): print_kpts_unblocked_upper(ints,fname,thresh_mono) - # ___ _ - # | ._ _|_ _ _ ._ _. | _ |_) o - # _|_ | | |_ (/_ (_| | (_| | _> |_) | - # _| - # + ########################################## + # # + # k-points # + # # + ########################################## + kconserv = tools.get_kconserv(cell, kpts) - qph5.create_dataset('nuclei/kconserv',data=np.transpose(kconserv+1,(0,2,1))) + + 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 - qph5.create_dataset('nuclei/kconserv_test',data=kcon_test) + with h5py.File(qph5path,'a') as qph5: + qph5.create_dataset('nuclei/kconserv_test',data=kcon_test) - - with open('K.qp','w') as outfile: - for a in range(Nk): - for b in range(Nk): - for c in range(Nk): - d = kconserv[a,b,c] - outfile.write('%s %s %s %s\n' % (a+1,c+1,b+1,d+1)) - + print_kcon_chem_to_phys(kconserv,'K.qp') + + ########################################## + # # + # Integrals Bi # + # # + ########################################## - intfile=h5py.File(mf.with_df._cderi,'r') +# qph5['ao_two_e_ints'].attrs['df_num']=naux - j3c = intfile.get('j3c') - naosq = nao*nao - naotri = (nao*(nao+1))//2 - j3ckeys = list(j3c.keys()) - j3ckeys.sort(key=lambda strkey:int(strkey)) - - # in new(?) version of PySCF, there is an extra layer of groups before the datasets - # datasets used to be [/j3c/0, /j3c/1, /j3c/2, ...] - # datasets now are [/j3c/0/0, /j3c/1/0, /j3c/2/0, ...] - j3clist = [j3c.get(i+'/0') for i in j3ckeys] - if j3clist==[None]*len(j3clist): - # if using older version, stop before last level - j3clist = [j3c.get(i) for i in j3ckeys] - - nkinvsq = 1./np.sqrt(Nk) - - # dimensions are (kikj,iaux,jao,kao), where kikj is compound index of kpts i and j - # output dimensions should be reversed (nao, nao, naux, nkptpairs) - j3arr=np.array([(i.value.reshape([-1,nao,nao]) if (i.shape[1] == naosq) else makesq3(i.value,nao)) * nkinvsq for i in j3clist]) + j3arr = get_j3ao(mf.with_df._cderi,nao,Nk) nkpt_pairs = j3arr.shape[0] naux = max(i.shape[0] for i in j3arr) print("n df fitting functions", naux) - qph5['ao_two_e_ints'].attrs['df_num']=naux - df_ao_tmp = np.zeros((nao,nao,naux,nkpt_pairs),dtype=np.complex128) + with h5py.File(qph5path,'a') as qph5: + qph5.create_group('ao_two_e_ints') + qph5['ao_two_e_ints'].attrs['df_num']=naux if print_ao_ints_df: - with open('D.qp','w') as outfile: - pass - with open('D.qp','a') as outfile: - for k,kpt_pair in enumerate(j3arr): - for iaux,dfbasfunc in enumerate(kpt_pair): - for i,i0 in enumerate(dfbasfunc): - for j,v in enumerate(i0): - if (abs(v) > bielec_int_threshold): - outfile.write(stri4z(i+1,j+1,iaux+1,k+1,v.real,v.imag)+'\n') - df_ao_tmp[i,j,iaux,k]=v + print_df(j3arr,'D.qp',bielec_int_threshold) + + df_ao_tmp = np.zeros((nao,nao,naux,nkpt_pairs),dtype=np.complex128) + for i,di in enumerate(j3arr): + df_ao_tmp[:,:,:di.shape[0],i] = np.transpose(di,(1,2,0)) - qph5.create_dataset('ao_two_e_ints/df_ao_integrals_real',data=df_ao_tmp.real) - qph5.create_dataset('ao_two_e_ints/df_ao_integrals_imag',data=df_ao_tmp.imag) + #df_ao_old = df_pad_ref_test(j3arr,nao,naux,nkpt_pairs) + #assert(abs(df_ao_tmp - df_ao_old).max() <= 1e-12) + + with h5py.File(qph5path,'a') as qph5: + qph5.create_dataset('ao_two_e_ints/df_ao_integrals_real',data=df_ao_tmp.real) + qph5.create_dataset('ao_two_e_ints/df_ao_integrals_imag',data=df_ao_tmp.imag) if print_mo_ints_df: - kpair_list=[] - for i in range(Nk): - for j in range(Nk): - if(i>=j): - kpair_list.append((i,j,idx2_tri((i,j)))) + from itertools import product + # WARNING: this is a generator, not a list; don't use it more than once + kpair_list = ((i,j,idx2_tri((i,j))) for (i,j) in product(range(Nk),repeat=2) if (i>=j)) j3mo = np.array([np.einsum('mij,ik,jl->mkl',j3arr[kij],mo_k[ki].conj(),mo_k[kj]) for ki,kj,kij in kpair_list]) + print_df(j3mo,'D.mo.qp',bielec_int_threshold) + df_mo_tmp = np.zeros((nmo,nmo,naux,nkpt_pairs),dtype=np.complex128) - with open('D_mo.qp','w') as outfile: - pass - with open('D_mo.qp','a') as outfile: - for k,kpt_pair in enumerate(j3mo): - for iaux,dfbasfunc in enumerate(kpt_pair): - for i,i0 in enumerate(dfbasfunc): - for j,v in enumerate(i0): - if (abs(v) > bielec_int_threshold): - outfile.write(stri4z(i+1,j+1,iaux+1,k+1,v.real,v.imag)+'\n') - df_mo_tmp[i,j,iaux,k]=v - qph5.create_dataset('mo_two_e_ints/df_mo_integrals_real',data=df_mo_tmp.real) - qph5.create_dataset('mo_two_e_ints/df_mo_integrals_imag',data=df_mo_tmp.imag) + for i,di in enumerate(j3mo): + df_mo_tmp[:,:,:di.shape[0],i] = np.transpose(di,(1,2,0)) + + #df_mo_old = df_pad_ref_test(j3mo,nmo,naux,nkpt_pairs) + #assert(abs(df_mo_tmp - df_mo_old).max() <= 1e-12) + + with h5py.File(qph5path,'a') as qph5: + qph5.create_dataset('mo_two_e_ints/df_mo_integrals_real',data=df_mo_tmp.real) + qph5.create_dataset('mo_two_e_ints/df_mo_integrals_imag',data=df_mo_tmp.imag) if (print_ao_ints_bi): print_ao_bi(mf,kconserv,'W.qp',bielec_int_threshold) @@ -835,7 +913,6 @@ def getj3ao(cell,mf, kpts, cas_idx=None, int_threshold = 1E-8): from pyscf.pbc import ao2mo from pyscf.pbc import tools from pyscf.pbc.gto import ecp - from pyscf.data import nist import h5py import scipy