From 02c6539daadd972760c06cd479008f749efe61c1 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 18 Feb 2020 14:11:22 -0600 Subject: [PATCH] fixed problem with iterating over unique 2-electron integrals should loop over union of two sets of integrals: set 1: i<=k j<=l ik<=jl set 2: i>k j kk)) cycle + if (ki>kl) cycle + ! if ((kl == kj) .and. (ki > kk)) cycle call idx2_tri_int(ki,kk,kikk2) - if (kikk2 > kjkl2) cycle + ! if (kikk2 > kjkl2) cycle if (ki < kk) then do i_ao=1,ao_num_per_kpt do j_ao=1,ao_num_per_kpt diff --git a/src/utils_complex/MolPyscfToQPkpts.py b/src/utils_complex/MolPyscfToQPkpts.py index e0284096..707599df 100644 --- a/src/utils_complex/MolPyscfToQPkpts.py +++ b/src/utils_complex/MolPyscfToQPkpts.py @@ -472,6 +472,199 @@ def pyscf2QP(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, outfile.write('%s %s %s %s %s %s\n' % (ii+1,jj+1,kk+1,ll+1,v.real,v.imag)) +def testprintbi(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8): + ''' + kpts = List of kpoints coordinates. Cannot be null, for gamma is other script + kmesh = Mesh of kpoints (optional) + cas_idx = List of active MOs. If not specified all MOs are actives + int_threshold = The integral will be not printed in they are bellow that + ''' + + 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 + + + bielec_int_threshold = int_threshold + + 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 + + + + # ___ _ + # | ._ _|_ _ _ ._ _. | _ |_) o + # _|_ | | |_ (/_ (_| | (_| | _> |_) | + # _| + # + kconserv = tools.get_kconserv(cell, kpts) + 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 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)) + + + intfile=h5py.File(mf.with_df._cderi,'r') + + 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]) + + nkpt_pairs = j3arr.shape[0] + df_ao_tmp = np.zeros((nao,nao,naux,nkpt_pairs),dtype=np.complex128) + + 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('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag)) + df_ao_tmp[i,j,iaux,k]=v + + 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)))) + 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]) + 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('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag)) + 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) + + + +# eri_4d_ao = np.zeros((Nk,nao,Nk,nao,Nk,nao,Nk,nao), dtype=np.complex) +# for d, kd in enumerate(kpts): +# for c, kc in enumerate(kpts): +# if c > d: break +# idx2_cd = idx2_tri(c,d) +# for b, kb in enumerate(kpts): +# if b > d: break +# a = kconserv[b,c,d] +# if idx2_tri(a,b) > idx2_cd: continue +# if ((c==d) and (a>b)): continue +# ka = kpts[a] +# v = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4) +# v *= 1./Nk +# eri_4d_ao[a,:,b,:,c,:,d] = v +# +# eri_4d_ao = eri_4d_ao.reshape([Nk*nao]*4) + + + with open('W.qp','w') as outfile: + pass + for d, kd in enumerate(kpts): + for c, kc in enumerate(kpts): + if c > d: break + idx2_cd = idx2_tri((c,d)) + for b, kb in enumerate(kpts): + if b > d: break + a = kconserv[b,c,d] + #if idx2_tri((a,b)) > idx2_cd: continue + if a>d: continue + #if ((c==d) and (a>b)): continue + ka = kpts[a] + + with open('W.qp','a') as outfile: + eri_4d_ao_kpt = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4) + eri_4d_ao_kpt *= 1./Nk + for l in range(nao): + ll=l+d*nao + for j in range(nao): + jj=j+c*nao + if jj>ll: break + idx2_jjll = idx2_tri((jj,ll)) + for k in range(nao): + kk=k+b*nao + if kk>ll: break + for i in range(nao): + ii=i+a*nao + if idx2_tri((ii,kk)) > idx2_jjll: break + if ((jj==ll) and (ii>kk)): break + v=eri_4d_ao_kpt[i,k,j,l] + if (abs(v) > bielec_int_threshold): + outfile.write('%s %s %s %s %s %s\n' % (ii+1,jj+1,kk+1,ll+1,v.real,v.imag)) + + + def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, print_ao_ints_bi=False, print_mo_ints_bi=False, @@ -754,8 +947,9 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, for b, kb in enumerate(kpts): if b > d: break a = kconserv[b,c,d] - if idx2_tri((a,b)) > idx2_cd: continue - if ((c==d) and (a>b)): continue + #if idx2_tri((a,b)) > idx2_cd: continue + if a > d: continue + #if ((c==d) and (a>b)): continue ka = kpts[a] if print_ao_ints_bi: diff --git a/src/utils_complex/dump_ao_2e_from_df.irp.f b/src/utils_complex/dump_ao_2e_from_df.irp.f index 2115a872..1ca00e09 100644 --- a/src/utils_complex/dump_ao_2e_from_df.irp.f +++ b/src/utils_complex/dump_ao_2e_from_df.irp.f @@ -52,9 +52,10 @@ subroutine run_ao_dump do kk=1,kl ki=kconserv(kl,kk,kj) - if ((kl == kj) .and. (ki > kk)) cycle + if (ki > kl) cycle + !if ((kl == kj) .and. (ki > kk)) cycle call idx2_tri_int(ki,kk,kikk2) - if (kikk2 > kjkl2) cycle + !if (kikk2 > kjkl2) cycle if (ki < kk) then do i_ao=1,ao_num_per_kpt do j_ao=1,ao_num_per_kpt @@ -72,7 +73,7 @@ subroutine run_ao_dump (1.d0,0.d0), ints_ik, ao_num_kpt_2, & ints_jl, ao_num_kpt_2, & (0.d0,0.d0), ints_ikjl, ao_num_kpt_2) - + print'((A),4(I4))','IJKL',ki,kj,kk,kl do il=1,ao_num_per_kpt l=il+(kl-1)*ao_num_per_kpt do ij=1,ao_num_per_kpt diff --git a/src/utils_complex/dump_ao_2e_from_df_all.irp.f b/src/utils_complex/dump_ao_2e_from_df_all.irp.f new file mode 100644 index 00000000..5ca67f11 --- /dev/null +++ b/src/utils_complex/dump_ao_2e_from_df_all.irp.f @@ -0,0 +1,117 @@ +program dump_ao_2e_from_df + call run_ao_dump +end + +subroutine run_ao_dump + use map_module + implicit none + BEGIN_DOC + ! fill ao bielec integral map using 3-index df integrals + END_DOC + + integer :: i,k,j,l + integer :: ki,kk,kj,kl + integer :: ii,ik,ij,il + integer :: kikk2,kjkl2,jl2,ik2 + integer :: i_ao,j_ao,i_df + + complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:), ints_ikjl(:,:,:,:) + + complex*16 :: integral,intmap, get_ao_two_e_integral_complex + double precision :: tmp_re,tmp_im + integer :: ao_num_kpt_2 + + logical :: use_map1 + integer(keY_kind) :: idx_tmp + double precision :: sign + + ao_num_kpt_2 = ao_num_per_kpt * ao_num_per_kpt + + + allocate( ints_jl(ao_num_per_kpt,ao_num_per_kpt,df_num)) + + do kl=1, kpt_num + do kj=1, kpt_num + call idx2_tri_int(kj,kl,kjkl2) + if (kj < kl) then + do i_ao=1,ao_num_per_kpt + do j_ao=1,ao_num_per_kpt + do i_df=1,df_num + ints_jl(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kjkl2)) + enddo + enddo + enddo + else + ints_jl = df_ao_integrals_complex(:,:,:,kjkl2) + endif + + allocate( & + ints_ik(ao_num_per_kpt,ao_num_per_kpt,df_num), & + ints_ikjl(ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt) & + ) + + do kk=1,kpt_num + ki=kconserv(kl,kk,kj) +! if ((kl == kj) .and. (ki > kk)) cycle + call idx2_tri_int(ki,kk,kikk2) +! if (kikk2 > kjkl2) cycle + if (ki < kk) then + do i_ao=1,ao_num_per_kpt + do j_ao=1,ao_num_per_kpt + do i_df=1,df_num + ints_ik(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kikk2)) + enddo + enddo + enddo +! ints_ik = conjg(reshape(df_mo_integral_array(:,:,:,kikk2),(/mo_num_per_kpt,mo_num_per_kpt,df_num/),order=(/2,1,3/))) + else + ints_ik = df_ao_integrals_complex(:,:,:,kikk2) + endif + + call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, df_num, & + (1.d0,0.d0), ints_ik, ao_num_kpt_2, & + ints_jl, ao_num_kpt_2, & + (0.d0,0.d0), ints_ikjl, ao_num_kpt_2) + print'((A),4(I4))','IJKL',ki,kj,kk,kl + do il=1,ao_num_per_kpt + l=il+(kl-1)*ao_num_per_kpt + do ij=1,ao_num_per_kpt + j=ij+(kj-1)*ao_num_per_kpt +! if (j>l) exit + call idx2_tri_int(j,l,jl2) + do ik=1,ao_num_per_kpt + k=ik+(kk-1)*ao_num_per_kpt +! if (k>l) exit + do ii=1,ao_num_per_kpt + i=ii+(ki-1)*ao_num_per_kpt +! if ((j==l) .and. (i>k)) exit +! call idx2_tri_int(i,k,ik2) +! if (ik2 > jl2) exit + integral = ints_ikjl(ii,ik,ij,il) + intmap = get_ao_two_e_integral_complex(i,j,k,l,ao_integrals_map,ao_integrals_map_2) +! print*,i,k,j,l,real(integral),imag(integral) + if ((cdabs(integral) + cdabs(intmap)) < ao_integrals_threshold) then + cycle + endif + if (cdabs(integral-intmap) < 1.d-14) then + cycle + !print'(4(I4),4(E15.7))',i,j,k,l,integral,intmap + else + print'(4(I4),4(E15.7),(A))',i,j,k,l,integral,intmap,'***' + endif + enddo !ii + enddo !ik + enddo !ij + enddo !il + enddo !kk + deallocate( & + ints_ik, & + ints_ikjl & + ) + enddo !kj + enddo !kl + deallocate( ints_jl ) + + +end +