diff --git a/.gitignore b/.gitignore index 096e385b..7b248722 100644 --- a/.gitignore +++ b/.gitignore @@ -3,7 +3,6 @@ quantum_package_static.tar.gz build.ninja .ninja_log .ninja_deps -bin/ lib/ config/qp_create_ninja.pickle src/*/.gitignore diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio index f980317b..f3e78d53 100755 --- a/bin/qp_convert_h5_to_ezfio +++ b/bin/qp_convert_h5_to_ezfio @@ -3,12 +3,13 @@ convert hdf5 output (e.g. from PySCF) to ezfio Usage: - qp_convert_h5_to_ezfio [--noqmc] [-o EZFIO_DIR] FILE + 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 @@ -20,6 +21,104 @@ 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 iilkjm", 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'] @@ -184,16 +882,16 @@ def convert_kpts(filename,qph5path,qmcpack=True): 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 + 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) @@ -202,9 +900,9 @@ def convert_kpts(filename,qph5path,qmcpack=True): 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 @@ -214,57 +912,57 @@ def convert_kpts(filename,qph5path,qmcpack=True): # #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(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 # @@ -302,7 +1000,7 @@ def convert_kpts(filename,qph5path,qmcpack=True): # 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) # @@ -319,7 +1017,7 @@ def convert_kpts(filename,qph5path,qmcpack=True): 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']) @@ -337,15 +1035,15 @@ def convert_kpts(filename,qph5path,qmcpack=True): 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() @@ -353,63 +1051,63 @@ def convert_kpts(filename,qph5path,qmcpack=True): 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(): @@ -428,7 +1126,7 @@ def convert_kpts(filename,qph5path,qmcpack=True): #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) @@ -439,13 +1137,13 @@ def convert_kpts(filename,qph5path,qmcpack=True): 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'] @@ -453,25 +1151,25 @@ def convert_cplx(filename,qph5path): 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 + 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 @@ -481,56 +1179,56 @@ def convert_cplx(filename,qph5path): # #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 # @@ -568,78 +1266,79 @@ def convert_cplx(filename,qph5path): 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(): @@ -653,7 +1352,7 @@ def convert_cplx(filename,qph5path): 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) @@ -664,15 +1363,22 @@ def convert_cplx(filename,qph5path): 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: @@ -680,11 +1386,20 @@ if __name__ == '__main__': if ARGUMENTS["--noqmc"]: qmcpack = False + if ARGUMENTS["--rmg"]: + rmg = True with h5py.File(FILE,'r') as qph5: - do_kpts = ('kconserv' in qph5['nuclei'].keys()) - if (do_kpts): + 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") - convert_kpts(EZFIO_FILE,FILE,qmcpack) + 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) @@ -697,4 +1412,3 @@ if __name__ == '__main__': # #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).""") - diff --git a/src/ao_two_e_ints/EZFIO.cfg b/src/ao_two_e_ints/EZFIO.cfg index 5b50f718..aab41e1f 100644 --- a/src/ao_two_e_ints/EZFIO.cfg +++ b/src/ao_two_e_ints/EZFIO.cfg @@ -35,3 +35,26 @@ doc: Real part of the df integrals over AOs size: (2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num) interface: ezfio +[chol_num] +type: integer +doc: number of cholesky vecs for each kpt +size: (nuclei.unique_kpt_num) +interface: ezfio, provider + +[chol_num_max] +type: integer +doc: max number of cholesky vecs +interface: ezfio, provider + +[io_chol_ao_integrals] +type: Disk_access +doc: Read/Write chol |AO| integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[chol_ao_integrals_complex] +type: double precision +doc: Cholesky decomposed integrals over AOs +size: (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) +interface: ezfio + diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f new file mode 100644 index 00000000..73fee7c7 --- /dev/null +++ b/src/ao_two_e_ints/cd_ao_ints.irp.f @@ -0,0 +1,261 @@ +!BEGIN_PROVIDER [ integer, chol_num_max ] +! implicit none +! BEGIN_DOC +! ! Max number of cholesky vectors. +! END_DOC +! chol_num_max = maxval(chol_num) +!END_PROVIDER + +BEGIN_PROVIDER [complex*16, chol_ao_integrals_complex, (ao_num_per_kpt,ao_num_per_kpt,chol_num_max,kpt_num,unique_kpt_num)] + implicit none + BEGIN_DOC + ! CD AO integrals + ! first two dims are AOs x AOs + ! 3rd dim is chol_vec (pad with zeros to max size to avoid dealing with ragged array) + ! 4th dim is over all kpts + ! last dim is over "unique" kpts (one for each pair of additive inverses modulo G) + END_DOC + integer :: i,j,k,l + + if (read_chol_ao_integrals) then + call ezfio_get_ao_two_e_ints_chol_ao_integrals_complex(chol_ao_integrals_complex) + print *, 'CD AO integrals read from disk' + else + print*,'CD AO integrals must be provided',irp_here + stop -1 + endif + + if (write_chol_ao_integrals) then + call ezfio_set_ao_two_e_ints_chol_ao_integrals_complex(chol_ao_integrals_complex) + print *, 'CD AO integrals written to disk' + endif + +END_PROVIDER + + +subroutine ao_map_fill_from_chol + use map_module + implicit none + BEGIN_DOC + ! TODO: check indexing/conj.transp. of slices; restructure loops + ! fill ao bielec integral map using 3-index cd 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_cd,kq + + complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:), ints_ikjl(:,:,:,:) + + complex*16 :: integral + integer :: n_integrals_1, n_integrals_2 + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i_1(:), buffer_i_2(:) + real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:) + double precision :: tmp_re,tmp_im + integer :: ao_num_kpt_2 + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + double precision :: map_mb + + logical :: use_map1 + integer(keY_kind) :: idx_tmp + double precision :: sign + + ao_num_kpt_2 = ao_num_per_kpt * ao_num_per_kpt + + size_buffer = min(ao_num_per_kpt*ao_num_per_kpt*ao_num_per_kpt,16000000) + print*, 'Providing the ao_bielec integrals from 3-index cholesky integrals' + call write_time(6) +! call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write') +! TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals + + call wall_time(wall_1) + call cpu_time(cpu_1) + + allocate( ints_jl(ao_num_per_kpt,ao_num_per_kpt,chol_num_max)) + + wall_0 = wall_1 + !TODO: change loops so that we only iterate over "correct" slices (i.e. ik block is stored directly, not as conj. transp.) + ! possible cases for (ik,jl) are (+,+), (+,-), (-,+), (-,-) + ! where + is the slice used as stored, and - is the conj. transp. of the stored data + ! (+,+) and (-,-) give the same information; we should always use (+,+) + ! (+,-) and (-,+) give the same information; we should always use (+,-) + do kQ = 1, kpt_num + do kl = 1, kpt_num + kj = qktok2(kQ,kl) + assert(kQ == qktok2(kj,kl)) + if (kj>kl) cycle + call idx2_tri_int(kj,kl,kjkl2) + !TODO: verify the kj, kl as 4th index in expressions below + if (kpt_sparse_map(kQ) > 0) then + !ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ)) + ints_jl = dconjg(chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ))) + else + do i_ao=1,ao_num_per_kpt + do j_ao=1,ao_num_per_kpt + do i_cd=1,chol_num_max + !ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ))) + ints_jl(i_ao,j_ao,i_cd) = chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ)) + enddo + enddo + enddo + endif + + !$OMP PARALLEL PRIVATE(i,k,j,l,ki,kk,ii,ik,ij,il,kikk2,jl2,ik2, & + !$OMP ints_ik, ints_ikjl, i_ao, j_ao, i_cd, & + !$OMP n_integrals_1, buffer_i_1, buffer_values_1, & + !$OMP n_integrals_2, buffer_i_2, buffer_values_2, & + !$OMP idx_tmp, tmp_re, tmp_im, integral,sign,use_map1) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer, kpt_num, ao_num_per_kpt, ao_num_kpt_2, & + !$OMP chol_num_max, chol_num, unique_kpt_num, kpt_sparse_map, qktok2, minusk, & + !$OMP kl,kj,kjkl2,ints_jl,kQ, & + !$OMP kconserv, chol_ao_integrals_complex, ao_integrals_threshold, ao_integrals_map, ao_integrals_map_2) + + allocate( & + ints_ik(ao_num_per_kpt,ao_num_per_kpt,chol_num_max), & + ints_ikjl(ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt), & + buffer_i_1(size_buffer), & + buffer_i_2(size_buffer), & + buffer_values_1(size_buffer), & + buffer_values_2(size_buffer) & + ) + + !$OMP DO SCHEDULE(guided) + do kk=1,kl + !print*,'debug' + !print*,kQ,kl,kj,kk + ki = qktok2(minusk(kk),kQ) + assert(ki == kconserv(kl,kk,kj)) + if (ki>kl) cycle + ! if ((kl == kj) .and. (ki > kk)) cycle + call idx2_tri_int(ki,kk,kikk2) + ! if (kikk2 > kjkl2) cycle + !TODO: check this! (ki, kk slice index and transpose/notranspose) + if (kpt_sparse_map(kQ) > 0) then + ints_ik = chol_ao_integrals_complex(:,:,:,ki,kpt_sparse_map(kQ)) + else + do i_ao=1,ao_num_per_kpt + do j_ao=1,ao_num_per_kpt + do i_cd=1,chol_num_max + ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kk,-kpt_sparse_map(kQ))) + enddo + enddo + enddo + endif + + call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, chol_num(kQ), & + (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) + + n_integrals_1=0 + n_integrals_2=0 + 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) +! print*,i,k,j,l,real(integral),imag(integral) + if (cdabs(integral) < ao_integrals_threshold) then + cycle + endif + call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign) + tmp_re = dble(integral) + tmp_im = dimag(integral) + if (use_map1) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp + buffer_values_1(n_integrals_1)=tmp_re + if (sign.ne.0.d0) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp+1 + buffer_values_1(n_integrals_1)=tmp_im*sign + endif + if (n_integrals_1 >= size(buffer_i_1)-1) then + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + n_integrals_1 = 0 + endif + else + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp + buffer_values_2(n_integrals_2)=tmp_re + if (sign.ne.0.d0) then + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp+1 + buffer_values_2(n_integrals_2)=tmp_im*sign + endif + if (n_integrals_2 >= size(buffer_i_2)-1) then + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + n_integrals_2 = 0 + endif + endif + + enddo !ii + enddo !ik + enddo !ij + enddo !il + + if (n_integrals_1 > 0) then + call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + endif + if (n_integrals_2 > 0) then + call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + endif + enddo !kk + !$OMP END DO NOWAIT + deallocate( & + ints_ik, & + ints_ikjl, & + buffer_i_1, & + buffer_i_2, & + buffer_values_1, & + buffer_values_2 & + ) + !$OMP END PARALLEL + enddo !kl + call wall_time(wall_2) + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(kQ)/float(kpt_num), '% in ', & + wall_2-wall_1,'s',map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + endif + + enddo !kQ + deallocate( ints_jl ) + + call map_sort(ao_integrals_map) + call map_unique(ao_integrals_map) + call map_sort(ao_integrals_map_2) + call map_unique(ao_integrals_map_2) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_1',ao_integrals_map) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_2',ao_integrals_map_2) + !call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + + call wall_time(wall_2) + call cpu_time(cpu_2) + + integer*8 :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() + + print*,'AO integrals provided:' + print*,' Size of AO map ', map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + print*,' Number of AO integrals: ', ao_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + +end subroutine ao_map_fill_from_chol + diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index c2a48a2f..662492b8 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -358,6 +358,11 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] print*, 'AO integrals provided from 3-index ao ints (periodic)' ao_two_e_integrals_in_map = .True. return + else if (read_chol_ao_integrals) then + call ao_map_fill_from_chol + print*, 'AO integrals provided from 3-index Cholesky ao ints (periodic)' + ao_two_e_integrals_in_map = .True. + return else print*,'calculation of periodic AOs not implemented' stop -1 diff --git a/src/hartree_fock/hf_energy.irp.f b/src/hartree_fock/hf_energy.irp.f index 5a68164f..ea5e4b1c 100644 --- a/src/hartree_fock/hf_energy.irp.f +++ b/src/hartree_fock/hf_energy.irp.f @@ -13,19 +13,22 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, hf_energy] &BEGIN_PROVIDER [ double precision, hf_two_electron_energy] +&BEGIN_PROVIDER [ double precision, hf_two_electron_energy_jk, (2)] &BEGIN_PROVIDER [ double precision, hf_one_electron_energy] implicit none BEGIN_DOC ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. END_DOC - integer :: i,j,k + integer :: i,j,k,jk hf_energy = nuclear_repulsion hf_two_electron_energy = 0.d0 + hf_two_electron_energy_jk = 0.d0 hf_one_electron_energy = 0.d0 if (is_complex) then - complex*16 :: hf_1e_tmp, hf_2e_tmp + complex*16 :: hf_1e_tmp, hf_2e_tmp, hf_2e_tmp_jk(2) hf_1e_tmp = (0.d0,0.d0) hf_2e_tmp = (0.d0,0.d0) + hf_2e_tmp_jk = (0.d0,0.d0) do k=1,kpt_num do j=1,ao_num_per_kpt do i=1,ao_num_per_kpt @@ -33,9 +36,21 @@ END_PROVIDER +ao_two_e_integral_beta_kpts(i,j,k) * scf_density_matrix_ao_beta_kpts(j,i,k) ) hf_1e_tmp += ao_one_e_integrals_kpts(i,j,k) * (scf_density_matrix_ao_alpha_kpts(j,i,k) & + scf_density_matrix_ao_beta_kpts (j,i,k) ) + do jk=1,2 + hf_2e_tmp_jk(jk) += 0.5d0 * ( ao_two_e_integral_alpha_kpts_jk(i,j,k,jk) * scf_density_matrix_ao_alpha_kpts(j,i,k) & + +ao_two_e_integral_beta_kpts_jk(i,j,k,jk) * scf_density_matrix_ao_beta_kpts(j,i,k) ) + enddo enddo enddo enddo + do jk=1,2 + if (dabs(dimag(hf_2e_tmp_jk(jk))).gt.1.d-10) then + print*,'HF_2e energy (jk) should be real:',jk,irp_here + stop -1 + else + hf_two_electron_energy_jk(jk) = dble(hf_2e_tmp_jk(jk)) + endif + enddo if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then print*,'HF_2e energy should be real:',irp_here stop -1 diff --git a/src/hartree_fock/print_e_scf.irp.f b/src/hartree_fock/print_e_scf.irp.f index 989c0b9c..8efeba0d 100644 --- a/src/hartree_fock/print_e_scf.irp.f +++ b/src/hartree_fock/print_e_scf.irp.f @@ -15,6 +15,9 @@ subroutine run print*,hf_one_electron_energy print*,hf_two_electron_energy print*,hf_energy + print*,'hf 2e J,K energy' + print*,hf_two_electron_energy_jk(1) + print*,hf_two_electron_energy_jk(2) end diff --git a/src/mo_two_e_ints/EZFIO.cfg b/src/mo_two_e_ints/EZFIO.cfg index c708792f..e8094c65 100644 --- a/src/mo_two_e_ints/EZFIO.cfg +++ b/src/mo_two_e_ints/EZFIO.cfg @@ -29,3 +29,15 @@ doc: Complex df integrals over MOs size: (2,mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,ao_two_e_ints.df_num,nuclei.kpt_pair_num) interface: ezfio +[io_chol_mo_integrals] +type: Disk_access +doc: Read/Write chol |MO| integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[chol_mo_integrals_complex] +type: double precision +doc: Cholesky decomposed integrals over MOs +size: (2,mo_basis.mo_num_per_kpt,mo_basis.mo_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num) +interface: ezfio + diff --git a/src/mo_two_e_ints/cd_mo_ints.irp.f b/src/mo_two_e_ints/cd_mo_ints.irp.f new file mode 100644 index 00000000..a832d510 --- /dev/null +++ b/src/mo_two_e_ints/cd_mo_ints.irp.f @@ -0,0 +1,325 @@ +BEGIN_PROVIDER [complex*16, chol_mo_integrals_complex, (mo_num_per_kpt,mo_num_per_kpt,chol_num_max,kpt_num,unique_kpt_num)] + implicit none + BEGIN_DOC + ! CD MO integrals + END_DOC + integer :: i,j,k,l + + if (read_chol_mo_integrals) then + call ezfio_get_mo_two_e_ints_chol_mo_integrals_complex(chol_mo_integrals_complex) + print *, 'CD MO integrals read from disk' + else + call chol_mo_from_chol_ao(chol_mo_integrals_complex,chol_ao_integrals_complex,mo_num_per_kpt,ao_num_per_kpt, & + chol_num_max,kpt_num,unique_kpt_num) + endif + + if (write_chol_mo_integrals) then + call ezfio_set_mo_two_e_ints_chol_mo_integrals_complex(chol_mo_integrals_complex) + print *, 'CD MO integrals written to disk' + endif + +END_PROVIDER + +subroutine mo_map_fill_from_chol_dot + use map_module + implicit none + BEGIN_DOC + ! TODO: verify correct indexing and conj.transp. + ! fill mo bielec integral map using 3-index cd integrals + END_DOC + + integer :: i,k,j,l,mu + integer :: ki,kk,kj,kl + integer :: ii,ik,ij,il + integer :: kikk2,kjkl2,jl2,ik2 + integer :: i_mo,j_mo,i_cd + integer :: kQ, Q_idx + + complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:) + + complex*16 :: integral,mjl,mik + integer :: n_integrals_1, n_integrals_2 + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i_1(:), buffer_i_2(:) + real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:) + double precision :: tmp_re,tmp_im + integer :: mo_num_kpt_2 + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + double precision :: map_mb + + logical :: use_map1 + integer(key_kind) :: idx_tmp + double precision :: sign + !complex*16, external :: zdotc + complex*16, external :: zdotu + + mo_num_kpt_2 = mo_num_per_kpt * mo_num_per_kpt + + size_buffer = min(mo_num_per_kpt*mo_num_per_kpt*mo_num_per_kpt,16000000) + print*, 'Providing the mo_bielec integrals from 3-index CD integrals' + call write_time(6) +! call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write') +! TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals + + call wall_time(wall_1) + call cpu_time(cpu_1) + + allocate( ints_jl(chol_num_max,mo_num_per_kpt,mo_num_per_kpt)) + allocate( ints_ik(chol_num_max,mo_num_per_kpt,mo_num_per_kpt)) + + wall_0 = wall_1 + do kQ = 1, kpt_num + Q_idx = kpt_sparse_map(kQ) + do kl = 1, kpt_num + kj = qktok2(kQ,kl) + assert(kQ == qktok2(kj,kl)) + if (kj>kl) cycle + call idx2_tri_int(kj,kl,kjkl2) + ints_jl = 0.d0 + if (Q_idx > 0) then + do i_mo=1,mo_num_per_kpt + do j_mo=1,mo_num_per_kpt + do i_cd=1,chol_num(kQ) + !ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx) + ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx)) + enddo + enddo + enddo + else + do i_mo=1,mo_num_per_kpt + do j_mo=1,mo_num_per_kpt + do i_cd=1,chol_num(kQ) + !ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx)) + ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx) + enddo + enddo + enddo + endif + + do kk=1,kl + ki = qktok2(minusk(kk),kQ) + assert(ki == kconserv(kl,kk,kj)) + if (ki>kl) cycle + call idx2_tri_int(ki,kk,kikk2) + ints_ik = 0.d0 + if (Q_idx > 0) then + do i_mo=1,mo_num_per_kpt + do j_mo=1,mo_num_per_kpt + do i_cd=1,chol_num(kQ) + ints_ik(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(i_mo,j_mo,i_cd,ki,Q_idx) + 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 + do i_mo=1,mo_num_per_kpt + do j_mo=1,mo_num_per_kpt + do i_cd=1,chol_num(kQ) + ints_ik(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(j_mo,i_mo,i_cd,kk,-Q_idx)) + enddo + enddo + enddo + endif + + !$OMP PARALLEL PRIVATE(i,k,j,l,ii,ik,ij,il,jl2,ik2, & + !$OMP mu, mik, mjl, & + !$OMP n_integrals_1, buffer_i_1, buffer_values_1, & + !$OMP n_integrals_2, buffer_i_2, buffer_values_2, & + !$OMP idx_tmp, tmp_re, tmp_im, integral,sign,use_map1) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer, kpt_num, mo_num_per_kpt, mo_num_kpt_2, & + !$OMP kl,kj,kjkl2,ints_jl, & + !$OMP ki,kk,kikk2,ints_ik, & + !$OMP kQ, Q_idx, chol_num, & + !$OMP kconserv, chol_mo_integrals_complex, mo_integrals_threshold, & + !$OMP mo_integrals_map, mo_integrals_map_2) + + allocate( & + buffer_i_1(size_buffer), & + buffer_i_2(size_buffer), & + buffer_values_1(size_buffer), & + buffer_values_2(size_buffer) & + ) + + n_integrals_1=0 + n_integrals_2=0 + !$OMP DO SCHEDULE(guided) + do il=1,mo_num_per_kpt + l=il+(kl-1)*mo_num_per_kpt + do ij=1,mo_num_per_kpt + j=ij+(kj-1)*mo_num_per_kpt + if (j>l) exit + call idx2_tri_int(j,l,jl2) + do ik=1,mo_num_per_kpt + k=ik+(kk-1)*mo_num_per_kpt + if (k>l) exit + do ii=1,mo_num_per_kpt + i=ii+(ki-1)*mo_num_per_kpt + if ((j==l) .and. (i>k)) exit + call idx2_tri_int(i,k,ik2) + if (ik2 > jl2) exit + !integral = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) + !integral = zdotu(chol_num(kQ),ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1) + integral = zdotu(chol_num(kQ),ints_jl(1,il,ij),1,ints_ik(1,ii,ik),1) +! print*,i,k,j,l,real(integral),imag(integral) + if (cdabs(integral) < mo_integrals_threshold) then + cycle + endif + call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign) + tmp_re = dble(integral) + tmp_im = dimag(integral) + if (use_map1) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp + buffer_values_1(n_integrals_1)=tmp_re + if (sign.ne.0.d0) then + n_integrals_1 += 1 + buffer_i_1(n_integrals_1)=idx_tmp+1 + buffer_values_1(n_integrals_1)=tmp_im*sign + endif + if (n_integrals_1 >= size(buffer_i_1)-1) then + call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1) + !call insert_into_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold) + n_integrals_1 = 0 + endif + else + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp + buffer_values_2(n_integrals_2)=tmp_re + if (sign.ne.0.d0) then + n_integrals_2 += 1 + buffer_i_2(n_integrals_2)=idx_tmp+1 + buffer_values_2(n_integrals_2)=tmp_im*sign + endif + if (n_integrals_2 >= size(buffer_i_2)-1) then + call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2) + !call insert_into_mo_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2,mo_integrals_threshold) + n_integrals_2 = 0 + endif + endif + + enddo !ii + enddo !ik + enddo !ij + enddo !il + !$OMP END DO NOWAIT + + if (n_integrals_1 > 0) then + call map_append(mo_integrals_map, buffer_i_1, buffer_values_1, n_integrals_1) + !call insert_into_mo_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1,mo_integrals_threshold) + endif + if (n_integrals_2 > 0) then + call map_append(mo_integrals_map_2, buffer_i_2, buffer_values_2, n_integrals_2) + !call insert_into_mo_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2,mo_integrals_threshold) + endif + deallocate( & + buffer_i_1, & + buffer_i_2, & + buffer_values_1, & + buffer_values_2 & + ) + !$OMP END PARALLEL + enddo !kk + enddo !kl + call wall_time(wall_2) + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(kQ)/float(kpt_num), '% in ', & + wall_2-wall_1,'s',map_mb(mo_integrals_map),'+',map_mb(mo_integrals_map_2),'MB' + endif + + enddo !kQ + deallocate( ints_jl,ints_ik ) + + call map_sort(mo_integrals_map) + call map_unique(mo_integrals_map) + call map_sort(mo_integrals_map_2) + call map_unique(mo_integrals_map_2) + !call map_merge(mo_integrals_map) + !call map_merge(mo_integrals_map_2) + + !!call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_complex_1',mo_integrals_map) + !!call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_complex_2',mo_integrals_map_2) + !!call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read') + + call wall_time(wall_2) + call cpu_time(cpu_2) + + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + + print*,'MO integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_map),'+',map_mb(mo_integrals_map_2),'MB' + print*,' Number of MO integrals: ', mo_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + +end subroutine mo_map_fill_from_chol_dot + + +subroutine chol_mo_from_chol_ao(cd_mo,cd_ao,n_mo,n_ao,n_cd,n_k,n_unique_k) + use map_module + implicit none + BEGIN_DOC + ! create 3-idx mo ints from 3-idx ao ints + END_DOC + integer,intent(in) :: n_mo,n_ao,n_cd,n_k,n_unique_k + complex*16,intent(out) :: cd_mo(n_mo,n_mo,n_cd,n_k,n_unique_k) + complex*16,intent(in) :: cd_ao(n_ao,n_ao,n_cd,n_k,n_unique_k) + integer :: ki,kk,mu,kQ,Q_idx + complex*16,allocatable :: coef_i(:,:), coef_k(:,:), ints_ik(:,:), ints_tmp(:,:) + double precision :: wall_1,wall_2,cpu_1,cpu_2 + + print*,'providing 3-index CD MO integrals from 3-index CD AO integrals' + + cd_mo = 0.d0 + + call wall_time(wall_1) + call cpu_time(cpu_1) + allocate( & + coef_i(n_ao,n_mo),& + coef_k(n_ao,n_mo),& + ints_ik(n_ao,n_ao),& + ints_tmp(n_mo,n_ao)& + ) + + do ki=1, kpt_num + coef_i = mo_coef_complex_kpts(:,:,ki) + do kk=1, kpt_num + coef_k = mo_coef_complex_kpts(:,:,kk) + kQ = qktok2(kk,ki) + Q_idx = kpt_sparse_map(kQ) + if (Q_idx < 0) cycle + + do mu=1, chol_num(kQ) + ints_ik = cd_ao(:,:,mu,ki,Q_idx) + call zgemm('C','N',n_mo,n_ao,n_ao, & + (1.d0,0.d0), coef_i, n_ao, & + ints_ik, n_ao, & + (0.d0,0.d0), ints_tmp, n_mo) + + call zgemm('N','N',n_mo,n_mo,n_ao, & + (1.d0,0.d0), ints_tmp, n_mo, & + coef_k, n_ao, & + (0.d0,0.d0), cd_mo(:,:,mu,ki,Q_idx), n_mo) + enddo + enddo + call wall_time(wall_2) + print*,100.*float(ki)/kpt_num, '% in ', & + wall_2-wall_1, 's' + enddo + + deallocate( & + coef_i, & + coef_k, & + ints_ik, & + ints_tmp & + ) + call wall_time(wall_2) + call cpu_time(cpu_2) + print*,' 3-idx CD MO provided' + print*,' cpu time:',cpu_2-cpu_1,'s' + print*,' wall time:',wall_2-wall_1,'s ( x ',(cpu_2-cpu_1)/(wall_2-wall_1),')' + +end subroutine chol_mo_from_chol_ao diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index e2447d6b..68af93bb 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -47,6 +47,12 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] !call mo_map_fill_from_df_single call mo_map_fill_from_df_dot return + else if (read_chol_mo_integrals.or.read_chol_ao_integrals) then + PROVIDE chol_mo_integrals_complex + !call mo_map_fill_from_chol + !call mo_map_fill_from_chol_single + call mo_map_fill_from_chol_dot + return else PROVIDE ao_two_e_integrals_in_map endif diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index b4599b72..79774ea5 100644 --- a/src/nuclei/EZFIO.cfg +++ b/src/nuclei/EZFIO.cfg @@ -60,3 +60,32 @@ type: integer doc: array containing information about k-point symmetry size: (nuclei.kpt_num,nuclei.kpt_num,nuclei.kpt_num) interface: ezfio + +[qktok2] +type: integer +doc: mapping from pairs of kpts to total per electron +size: (nuclei.kpt_num,nuclei.kpt_num) +interface: ezfio + +[minusk] +type: integer +doc: additive inverse for each kpt +size: (nuclei.kpt_num) +interface: ezfio + +[kpt_sparse_map] +type: integer +doc: mapping from kpt idx to unique idx, negative for conj. transp. +size: (nuclei.kpt_num) +interface: ezfio + +[unique_kpt_num] +type: integer +doc: number of pairs of kpts that are additive inverses (mod G) +interface: ezfio, provider + +[io_kpt_symm] +doc: Read/Write kpt_symm arrays from/to disk [ Write | Read | None ] +type: Disk_access +interface: ezfio,provider,ocaml +default: None diff --git a/src/nuclei/kconserv_cplx.irp.f b/src/nuclei/kconserv_cplx.irp.f index 616ba779..e0958bec 100644 --- a/src/nuclei/kconserv_cplx.irp.f +++ b/src/nuclei/kconserv_cplx.irp.f @@ -21,8 +21,9 @@ BEGIN_PROVIDER [integer, kconserv, (kpt_num,kpt_num,kpt_num)] call ezfio_get_nuclei_kconserv(kconserv) print *, 'kconserv read from disk' else - print*,'kconserv must be provided' - stop -1 + call set_kconserv(kconserv) + !print*,'kconserv must be provided' + !stop -1 endif if (write_kconserv) then call ezfio_set_nuclei_kconserv(kconserv) @@ -30,6 +31,86 @@ BEGIN_PROVIDER [integer, kconserv, (kpt_num,kpt_num,kpt_num)] endif END_PROVIDER +BEGIN_PROVIDER [integer, qktok2, (kpt_num,kpt_num)] + implicit none + BEGIN_DOC + ! Information about k-point symmetry + ! + ! for k-points I,K: qktok2(K,I) = \alpha + ! where Q_{\alpha} = k_I - k_K + ! + END_DOC + + if (read_kpt_symm) then + call ezfio_get_nuclei_qktok2(qktok2) + print *, 'qktok2 read from disk' + else + print*,'qktok2 must be provided' + stop -1 + endif + if (write_kpt_symm) then + call ezfio_set_nuclei_qktok2(qktok2) + print *, 'qktok2 written to disk' + endif +END_PROVIDER + +BEGIN_PROVIDER [integer, minusk, (kpt_num)] + implicit none + BEGIN_DOC + ! Information about k-point symmetry + ! + ! for k-point I: minusk(I) = K + ! where k_I + k_K = 0 (mod G) + ! + END_DOC + + if (read_kpt_symm) then + call ezfio_get_nuclei_minusk(minusk) + print *, 'minusk read from disk' + else + print*,'minusk must be provided' + stop -1 + endif + if (write_kpt_symm) then + call ezfio_set_nuclei_minusk(minusk) + print *, 'minusk written to disk' + endif +END_PROVIDER + +BEGIN_PROVIDER [integer, kpt_sparse_map, (kpt_num)] + implicit none + BEGIN_DOC + ! Information about k-point symmetry + ! + ! for k-point I: if kpt_sparse_map(I) = j + ! if j>0: data for k_I is stored at index j in chol_ints + ! if j<0: data for k_I is conj. transp. of data at index j in chol_{ao,mo}_integrals_complex + ! + ! if we have h5 data stored under L[i]: + ! count=1 + ! do i=1,N_L + ! kpt_sparse_map(i)=count + ! if (minusk(i) != i) then + ! kpt_sparse_map(minusk(i)) = -count + ! endif + ! count += 1 + ! enddo + ! + END_DOC + + if (read_kpt_symm) then + call ezfio_get_nuclei_kpt_sparse_map(kpt_sparse_map) + print *, 'kpt_sparse_map read from disk' + else + print*,'kpt_sparse_map must be provided' + stop -1 + endif + if (write_kpt_symm) then + call ezfio_set_nuclei_kpt_sparse_map(kpt_sparse_map) + print *, 'kpt_sparse_map written to disk' + endif +END_PROVIDER + subroutine double_allowed_kpts(kh1,kh2,kp1,kp2,is_allowed) implicit none integer, intent(in) :: kh1,kh2,kp1,kp2 @@ -38,3 +119,19 @@ subroutine double_allowed_kpts(kh1,kh2,kp1,kp2,is_allowed) is_allowed = (kconserv(kh1,kh2,kp1) == kp2) end subroutine +subroutine set_kconserv(kcon) + implicit none + integer, intent(out) :: kcon(kpt_num,kpt_num,kpt_num) + integer :: i,j,k,qik + + do i=1,kpt_num + do k=1,kpt_num + ! Q = k_I - k_K + qik = qktok2(k,i) + do j=1,kpt_num + ! k_L = k_J - (-(k_I - k_K)) + kcon(i,j,k) = qktok2(minusk(j),qik) + enddo + enddo + enddo +end subroutine diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index e2ada6fc..2f7d4a1a 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/src/scf_utils/fock_matrix_cplx.irp.f @@ -542,27 +542,26 @@ BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_kpts, (ao_num_per_kpt, ao_num_per_kp endif END_PROVIDER - - BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] -&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts , (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts_jk, (ao_num_per_kpt, ao_num_per_kpt, kpt_num, 2) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts_jk , (ao_num_per_kpt, ao_num_per_kpt, kpt_num, 2) ] use map_module implicit none BEGIN_DOC - ! Alpha and Beta Fock matrices in AO basis set + ! Alpha and Beta Fock matrices in AO basis set separated into j/k END_DOC !TODO: finish implementing this: see complex qp1 (different mapping) - + integer :: i,j,k,l,k1,r,s integer :: i0,j0,k0,l0 integer*8 :: p,q complex*16 :: integral, c0 - complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:) - complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:) - - ao_two_e_integral_alpha_kpts = (0.d0,0.d0) - ao_two_e_integral_beta_kpts = (0.d0,0.d0) + complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:,:) + complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:,:) + + ao_two_e_integral_alpha_kpts_jk = (0.d0,0.d0) + ao_two_e_integral_beta_kpts_jk = (0.d0,0.d0) PROVIDE ao_two_e_integrals_in_map scf_density_matrix_ao_alpha_kpts scf_density_matrix_ao_beta_kpts - + integer(omp_lock_kind) :: lck(ao_num) integer(map_size_kind) :: i8 integer :: ii(4), jj(4), kk(4), ll(4), k2 @@ -572,7 +571,267 @@ END_PROVIDER complex*16, parameter :: i_sign(4) = (/(0.d0,1.d0),(0.d0,1.d0),(0.d0,-1.d0),(0.d0,-1.d0)/) integer(key_kind) :: key1 integer :: kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l - + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & + !$OMP kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l, & + !$OMP c0,key1)& + !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts, kpt_num, irp_here, & + !$OMP SCF_density_matrix_ao_beta_kpts, & + !$OMP ao_integrals_map, ao_two_e_integral_alpha_kpts_jk, ao_two_e_integral_beta_kpts_jk) + + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2), & + ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2)) + ao_two_e_integral_alpha_tmp = (0.d0,0.d0) + ao_two_e_integral_beta_tmp = (0.d0,0.d0) + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) + do k1=1,n_elements + ! get original key + ! reverse of 2*key (imag part) and 2*key-1 (real part) + key1 = shiftr(keys(k1)+1,1) + + call two_e_integrals_index_reverse_complex_1(ii,jj,kk,ll,key1) + ! i<=k, j<=l, ik<=jl + ! ijkl, jilk, klij*, lkji* + + if (shiftl(key1,1)==keys(k1)) then !imaginary part (even) + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate + + !G_a(i,k) += D_{ab}(l,j)*() + !G_b(i,k) += D_{ab}(l,j)*() + !G_a(i,l) -= D_a (k,j)*() + !G_b(i,l) -= D_b (k,j)*() + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + else ! real part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = values(k1) + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_kpts_jk += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_kpts_jk += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & + !$OMP kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l, & + !$OMP c0,key1)& + !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts,kpt_num, irp_here, & + !$OMP SCF_density_matrix_ao_beta_kpts, & + !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_kpts_jk, ao_two_e_integral_beta_kpts_jk) + + call get_cache_map_n_elements_max(ao_integrals_map_2,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2), & + ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2)) + ao_two_e_integral_alpha_tmp = (0.d0,0.d0) + ao_two_e_integral_beta_tmp = (0.d0,0.d0) + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map_2%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map_2,i8,keys,values,n_elements) + do k1=1,n_elements + ! get original key + ! reverse of 2*key (imag part) and 2*key-1 (real part) + key1 = shiftr(keys(k1)+1,1) + + call two_e_integrals_index_reverse_complex_2(ii,jj,kk,ll,key1) + ! i>=k, j<=l, ik<=jl + ! ijkl, jilk, klij*, lkji* + if (shiftl(key1,1)==keys(k1)) then !imaginary part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate + + !G_a(i,k) += D_{ab}(l,j)*() + !G_b(i,k) += D_{ab}(l,j)*() + !G_a(i,l) -= D_a (k,j)*() + !G_b(i,l) -= D_b (k,j)*() + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + else ! real part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = values(k1) + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_kpts_jk += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_kpts_jk += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + +END_PROVIDER + + + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts , (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] + use map_module + implicit none + BEGIN_DOC + ! Alpha and Beta Fock matrices in AO basis set + END_DOC + !TODO: finish implementing this: see complex qp1 (different mapping) + + integer :: i,j,k,l,k1,r,s + integer :: i0,j0,k0,l0 + integer*8 :: p,q + complex*16 :: integral, c0 + complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:) + complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:) + + ao_two_e_integral_alpha_kpts = (0.d0,0.d0) + ao_two_e_integral_beta_kpts = (0.d0,0.d0) + PROVIDE ao_two_e_integrals_in_map scf_density_matrix_ao_alpha_kpts scf_density_matrix_ao_beta_kpts + + integer(omp_lock_kind) :: lck(ao_num) + integer(map_size_kind) :: i8 + integer :: ii(4), jj(4), kk(4), ll(4), k2 + integer(cache_map_size_kind) :: n_elements_max, n_elements + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + complex*16, parameter :: i_sign(4) = (/(0.d0,1.d0),(0.d0,1.d0),(0.d0,-1.d0),(0.d0,-1.d0)/) + integer(key_kind) :: key1 + integer :: kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & @@ -581,14 +840,14 @@ END_PROVIDER !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts, kpt_num, irp_here, & !$OMP SCF_density_matrix_ao_beta_kpts, & !$OMP ao_integrals_map, ao_two_e_integral_alpha_kpts, ao_two_e_integral_beta_kpts) - + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) allocate(keys(n_elements_max), values(n_elements_max)) allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num), & ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num)) ao_two_e_integral_alpha_tmp = (0.d0,0.d0) ao_two_e_integral_beta_tmp = (0.d0,0.d0) - + !$OMP DO SCHEDULE(static,1) do i8=0_8,ao_integrals_map%map_size n_elements = n_elements_max @@ -686,7 +945,7 @@ END_PROVIDER deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) !$OMP END PARALLEL - + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & @@ -695,14 +954,14 @@ END_PROVIDER !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts,kpt_num, irp_here, & !$OMP SCF_density_matrix_ao_beta_kpts, & !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_kpts, ao_two_e_integral_beta_kpts) - + call get_cache_map_n_elements_max(ao_integrals_map_2,n_elements_max) allocate(keys(n_elements_max), values(n_elements_max)) allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num), & ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num)) ao_two_e_integral_alpha_tmp = (0.d0,0.d0) ao_two_e_integral_beta_tmp = (0.d0,0.d0) - + !$OMP DO SCHEDULE(static,1) do i8=0_8,ao_integrals_map_2%map_size n_elements = n_elements_max diff --git a/src/utils_complex/dump_ao_2e_cplx.irp.f b/src/utils_complex/dump_ao_2e_cplx.irp.f index 2db5f614..ea36b876 100644 --- a/src/utils_complex/dump_ao_2e_cplx.irp.f +++ b/src/utils_complex/dump_ao_2e_cplx.irp.f @@ -15,19 +15,21 @@ subroutine run do k=1,ao_num do l=1,ao_num tmp_cmplx = get_ao_two_e_integral_complex(i,j,k,l,ao_integrals_map,ao_integrals_map_2) - print'(4(I4),2(E23.15))',i,j,k,l,tmp_cmplx + if (cdabs(tmp_cmplx) .gt. 1E-10) then + print'(4(I4),2(E23.15))',i,k,j,l,tmp_cmplx + endif enddo enddo enddo enddo - print*,'map1' - do i=0,ao_integrals_map%map_size - print*,i,ao_integrals_map%map(i)%value(:) - print*,i,ao_integrals_map%map(i)%key(:) - enddo - print*,'map2' - do i=0,ao_integrals_map_2%map_size - print*,i,ao_integrals_map_2%map(i)%value(:) - print*,i,ao_integrals_map_2%map(i)%key(:) - enddo + !print*,'map1' + !do i=0,ao_integrals_map%map_size + ! print*,i,ao_integrals_map%map(i)%value(:) + ! print*,i,ao_integrals_map%map(i)%key(:) + !enddo + !print*,'map2' + !do i=0,ao_integrals_map_2%map_size + ! print*,i,ao_integrals_map_2%map(i)%value(:) + ! print*,i,ao_integrals_map_2%map(i)%key(:) + !enddo end diff --git a/src/utils_complex/dump_cd_ints.irp.f b/src/utils_complex/dump_cd_ints.irp.f new file mode 100644 index 00000000..40adb2b4 --- /dev/null +++ b/src/utils_complex/dump_cd_ints.irp.f @@ -0,0 +1,29 @@ +program dump_cd_ksym + call run +end + +subroutine run + use map_module + implicit none + + integer ::q,k,n,i,j + double precision :: vr, vi + complex*16 :: v + print*,"chol_ao_integrals_complex q,k,n,i,j" + provide chol_ao_integrals_complex + do q = 1, unique_kpt_num + do k = 1, kpt_num + do n = 1, chol_num_max + do i = 1, ao_num_per_kpt + do j = 1, ao_num_per_kpt + v = chol_ao_integrals_complex(i,j,n,k,q) + vr = dble(v) + vi = dimag(v) + print '(5(I6,X),2(E25.15,X))', q, k, n, i, j, vr, vi + enddo + enddo + enddo + enddo + enddo + +end diff --git a/src/utils_complex/dump_cd_ksym.irp.f b/src/utils_complex/dump_cd_ksym.irp.f new file mode 100644 index 00000000..d5f9cfdb --- /dev/null +++ b/src/utils_complex/dump_cd_ksym.irp.f @@ -0,0 +1,47 @@ +program dump_cd_ksym + call run +end + +subroutine run + use map_module + implicit none + + integer ::i,j,k,l + integer(key_kind) :: idx + logical :: use_map1 + double precision :: sign + do i=1,5 + do j=1,5 + do k=1,5 + do l=1,5 + call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx,sign) + print'(4(I4,X),(L6),(I8),(F10.1))',i,j,k,l,use_map1,idx,sign + enddo + enddo + enddo + enddo + + provide qktok2 minusk kconserv + print*,'minusk' + do i=1,kpt_num + j = minusk(i) + print'(2(I4))',i,j + enddo + print*,'qktok2' + do i=1,kpt_num + do j=1,kpt_num + k = qktok2(i,j) + print'(3(I4))',i,j,k + enddo + enddo + print*,'kconserv' + do i=1,kpt_num + do j=1,kpt_num + do k=1,kpt_num + l = kconserv(i,j,k) + print'(4(I4))',i,j,k,l + enddo + enddo + enddo + +end diff --git a/src/utils_complex/test_cd_ksym.irp.f b/src/utils_complex/test_cd_ksym.irp.f new file mode 100644 index 00000000..be1433be --- /dev/null +++ b/src/utils_complex/test_cd_ksym.irp.f @@ -0,0 +1,234 @@ +program test_cd_ksym + call run +end + +subroutine run + use map_module + implicit none + + !integer ::i,j,k,l + + provide qktok2 minusk kconserv + !print*,'minusk' + !do i=1,kpt_num + ! j = minusk(i) + ! print'(2(I4))',i,j + !enddo + !print*,'qktok2' + !do i=1,kpt_num + ! do j=1,kpt_num + ! k = qktok2(i,j) + ! print'(3(I4))',i,j,k + ! enddo + !enddo + !print*,'kconserv' + !do i=1,kpt_num + ! do j=1,kpt_num + ! do k=1,kpt_num + ! l = kconserv(i,j,k) + ! print'(4(I4))',i,j,k,l + ! enddo + ! enddo + !enddo + + 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_cd,kq + + complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:), ints_ikjl(:,:,:,:) + + complex*16 :: integral + integer :: n_integrals_1, n_integrals_2 + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i_1(:), buffer_i_2(:) + real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:) + double precision :: tmp_re,tmp_im + integer :: ao_num_kpt_2 + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + double precision :: map_mb + + logical :: use_map1 + integer(keY_kind) :: idx_tmp + double precision :: sign + + ao_num_kpt_2 = ao_num_per_kpt * ao_num_per_kpt + + size_buffer = min(ao_num_per_kpt*ao_num_per_kpt*ao_num_per_kpt,16000000) + print*, 'Providing the ao_bielec integrals from 3-index cholesky integrals' + call write_time(6) +! call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write') +! TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals + + call wall_time(wall_1) + call cpu_time(cpu_1) + + !allocate( ints_jl(ao_num_per_kpt,ao_num_per_kpt,chol_num_max)) + + wall_0 = wall_1 + ! ki + kj == kk + kl required for to be nonzero + !TODO: change loops so that we only iterate over "correct" slices (i.e. ik block is stored directly, not as conj. transp.) + ! possible cases for (ik,jl) are (+,+), (+,-), (-,+), (-,-) + ! where + is the slice used as stored, and - is the conj. transp. of the stored data + ! (+,+) and (-,-) give the same information; we should always use (+,+) + ! (+,-) and (-,+) give the same information; we should always use (+,-) + do kQ = 1, kpt_num + do kl = 1, kpt_num + kj = qktok2(kQ,kl) + assert(kQ == qktok2(kj,kl)) + if (kj>kl) cycle + call idx2_tri_int(kj,kl,kjkl2) + !TODO: verify the kj, kl as 4th index in expressions below + !if (kpt_sparse_map(kQ) > 0) then + ! ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ)) + !else + ! !do i_ao=1,ao_num_per_kpt + ! ! do j_ao=1,ao_num_per_kpt + ! ! do i_cd=1,chol_num_max + ! ! ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ))) + ! ! enddo + ! ! enddo + ! !enddo + !endif + + !allocate( & + ! ints_ik(ao_num_per_kpt,ao_num_per_kpt,chol_num_max), & + ! ints_ikjl(ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt), & + ! buffer_i_1(size_buffer), & + ! buffer_i_2(size_buffer), & + ! buffer_values_1(size_buffer), & + ! buffer_values_2(size_buffer) & + !) + + do kk=1,kl + ki = qktok2(minusk(kk),kQ) + assert(ki == kconserv(kl,kk,kj)) + if (ki>kl) cycle + ! if ((kl == kj) .and. (ki > kk)) cycle + call idx2_tri_int(ki,kk,kikk2) + print*,kQ,kl,kj,kk,ki + ! if (kikk2 > kjkl2) cycle + !TODO: check this! (ki, kk slice index and transpose/notranspose) + !if (kpt_sparse_map(kQ) > 0) then + ! ints_ik = chol_ao_integrals_complex(:,:,:,ki,kpt_sparse_map(kQ)) + !else + ! do i_ao=1,ao_num_per_kpt + ! do j_ao=1,ao_num_per_kpt + ! do i_cd=1,chol_num_max + ! ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kk,-kpt_sparse_map(kQ))) + ! enddo + ! enddo + ! enddo + !endif + + !call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, chol_num(kQ), & + ! (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) + + !n_integrals_1=0 + !n_integrals_2=0 + !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) +! ! print*,i,k,j,l,real(integral),imag(integral) + ! if (cdabs(integral) < ao_integrals_threshold) then + ! cycle + ! endif + ! call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign) + ! tmp_re = dble(integral) + ! tmp_im = dimag(integral) + ! !if (use_map1) then + ! ! n_integrals_1 += 1 + ! ! buffer_i_1(n_integrals_1)=idx_tmp + ! ! buffer_values_1(n_integrals_1)=tmp_re + ! ! if (sign.ne.0.d0) then + ! ! n_integrals_1 += 1 + ! ! buffer_i_1(n_integrals_1)=idx_tmp+1 + ! ! buffer_values_1(n_integrals_1)=tmp_im*sign + ! ! endif + ! ! if (n_integrals_1 >= size(buffer_i_1)-1) then + ! ! call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + ! ! n_integrals_1 = 0 + ! ! endif + ! !else + ! !n_integrals_2 += 1 + ! !buffer_i_2(n_integrals_2)=idx_tmp + ! !buffer_values_2(n_integrals_2)=tmp_re + ! !if (sign.ne.0.d0) then + ! ! n_integrals_2 += 1 + ! ! buffer_i_2(n_integrals_2)=idx_tmp+1 + ! ! buffer_values_2(n_integrals_2)=tmp_im*sign + ! !endif + ! !if (n_integrals_2 >= size(buffer_i_2)-1) then + ! ! call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + ! ! n_integrals_2 = 0 + ! !endif + ! endif + + ! enddo !ii + ! enddo !ik + ! enddo !ij + !enddo !il + + !if (n_integrals_1 > 0) then + ! call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1) + !endif + !if (n_integrals_2 > 0) then + ! call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2) + !endif + enddo !kk + !deallocate( & + ! ints_ik, & + ! ints_ikjl, & + ! buffer_i_1, & + ! buffer_i_2, & + ! buffer_values_1, & + ! buffer_values_2 & + ! ) + enddo !kl + call wall_time(wall_2) + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + !print*, 100.*float(kQ)/float(kpt_num), '% in ', & + ! wall_2-wall_1,'s',map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + endif + + enddo !kQ + !deallocate( ints_jl ) + + !call map_sort(ao_integrals_map) + !call map_unique(ao_integrals_map) + !call map_sort(ao_integrals_map_2) + !call map_unique(ao_integrals_map_2) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_1',ao_integrals_map) + !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_2',ao_integrals_map_2) + !call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + + call wall_time(wall_2) + call cpu_time(cpu_2) + + !integer*8 :: get_ao_map_size, ao_map_size + !ao_map_size = get_ao_map_size() + + print*,'AO integrals provided:' + !print*,' Size of AO map ', map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB' + !print*,' Number of AO integrals: ', ao_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + +end