mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-19 04:22:32 +01:00
debugging cd ints
This commit is contained in:
Kevin Gasperich
parent
2371798bdc
commit
fe9ddc4d98
@ -21,7 +21,93 @@ from docopt import docopt
|
||||
import gzip
|
||||
#fname = sys.argv[1]
|
||||
#qph5name = sys.argv[2]
|
||||
def idx2_tri(i,j):
|
||||
"""
|
||||
for 0-indexed counting
|
||||
"""
|
||||
p = max(i,j)
|
||||
q = min(i,j)
|
||||
return q + (p*(p+1))//2
|
||||
def idx2_tri_1(i,j):
|
||||
"""
|
||||
for 1-indexed counting
|
||||
"""
|
||||
p = max(i,j)
|
||||
q = min(i,j)
|
||||
return q + (p*(p-1))//2
|
||||
|
||||
def idx4_cplx_1(i,j,k,l):
|
||||
"""
|
||||
original function from qp2 (fortran counting)
|
||||
"""
|
||||
p = idx2_tri_1(i,k)
|
||||
q = idx2_tri_1(j,l)
|
||||
i1 = idx2_tri_1(p,q)
|
||||
return (i1,p,q)
|
||||
|
||||
def ao_idx_map_sign(i,j,k,l):
|
||||
"""
|
||||
qp2 indexing
|
||||
"""
|
||||
idx,ik,jl = idx4_cplx_1(i,j,k,l)
|
||||
ij = idx2_tri_1(i,j)
|
||||
kl = idx2_tri_1(k,l)
|
||||
idx = 2*idx - 1
|
||||
|
||||
if (ij==kl):
|
||||
sign = 0.0
|
||||
use_map1 = False
|
||||
else:
|
||||
if ik==jl:
|
||||
if i<k:
|
||||
sign = 1.0
|
||||
use_map1 = True
|
||||
else:
|
||||
sign = -1.0
|
||||
use_map1 = True
|
||||
elif i==k:
|
||||
if j<l:
|
||||
sign = 1.0
|
||||
use_map1 = True
|
||||
else:
|
||||
sign = -1.0
|
||||
use_map1 = True
|
||||
elif j==l:
|
||||
if i<k:
|
||||
sign = 1.0
|
||||
use_map1 = True
|
||||
else:
|
||||
sign = -1.0
|
||||
use_map1 = True
|
||||
elif ((i<k) == (j<l)):
|
||||
if i<k:
|
||||
sign = 1.0
|
||||
use_map1 = True
|
||||
else:
|
||||
sign = -1.0
|
||||
use_map1 = True
|
||||
else:
|
||||
if ((j<l) == (ik<jl)):
|
||||
sign = 1.0
|
||||
use_map1 = False
|
||||
else:
|
||||
sign = -1.0
|
||||
use_map1 = False
|
||||
return (idx, use_map1, sign)
|
||||
def get_ao_int_cplx(i,j,k,l,map1,map2):
|
||||
idx,use_m1,sgn = ao_idx_map_sign(i,j,k,l)
|
||||
if use_m1:
|
||||
tmp_re = map1[idx]
|
||||
tmp_im = map1[idx+1]
|
||||
tmp_im *= sgn
|
||||
else:
|
||||
tmp_re = map2[idx]
|
||||
if sgn != 0.0:
|
||||
tmp_im = map2[idx+1]
|
||||
tmp_im *= sgn
|
||||
else:
|
||||
tmp_im = 0.0
|
||||
return tmp_re + 1j*tmp_im
|
||||
def kconserv_p_from_qkk2_mk(qkk2,mk):
|
||||
nk, nk2 = qkk2.shape
|
||||
assert(nk == nk2)
|
||||
@ -30,6 +116,7 @@ def kconserv_p_from_qkk2_mk(qkk2,mk):
|
||||
for j in range(nk):
|
||||
for k in range(nk):
|
||||
kcon_p[i,j,k] = qkk2[mk[j],qkk2[k,i]]
|
||||
assert(qkk2[mk[j],qkk2[k,i]] == qkk2[qkk2[j,k],i])
|
||||
return kcon_p
|
||||
|
||||
|
||||
@ -201,6 +288,11 @@ def convert_mol(filename,qph5path):
|
||||
def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
|
||||
import json
|
||||
from scipy.linalg import block_diag
|
||||
dump_fci, dump_cd = (False, False)
|
||||
#dump_fci, dump_cd = (True, False)
|
||||
#dump_fci, dump_cd = (False, True)
|
||||
#dump_fci2 = True
|
||||
dump_fci2 = False
|
||||
|
||||
ezfio.set_file(filename)
|
||||
ezfio.set_nuclei_is_complex(True)
|
||||
@ -385,10 +477,10 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
|
||||
isparse += 1
|
||||
else:
|
||||
kpt_sparse_map[i] = -kpt_sparse_map[minusk[i]-1]
|
||||
ezfio.set_nuclei_kconserv(kconserv)
|
||||
ezfio.set_nuclei_kconserv(kconserv.transpose(2,1,0))
|
||||
ezfio.set_nuclei_io_kconserv('Read')
|
||||
ezfio.set_nuclei_minusk(minusk)
|
||||
ezfio.set_nuclei_qktok2(QKTok2)
|
||||
ezfio.set_nuclei_qktok2(np.transpose(QKTok2)) #transposed for correct col-major ordering
|
||||
ezfio.set_nuclei_kpt_sparse_map(kpt_sparse_map)
|
||||
ezfio.set_nuclei_unique_kpt_num(unique_kpt_num)
|
||||
# kpt_sparse_map
|
||||
@ -445,7 +537,263 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
|
||||
#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)
|
||||
#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
|
||||
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}')
|
||||
|
||||
|
||||
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#)
|
||||
"""
|
||||
@ -991,6 +1339,12 @@ def convert_cplx(filename,qph5path):
|
||||
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
|
||||
|
||||
@ -92,11 +92,13 @@ subroutine ao_map_fill_from_chol
|
||||
!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
|
||||
|
||||
@ -81,7 +81,8 @@ subroutine mo_map_fill_from_chol_dot
|
||||
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) = 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
|
||||
@ -89,7 +90,8 @@ subroutine mo_map_fill_from_chol_dot
|
||||
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) = 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
|
||||
@ -268,7 +270,7 @@ subroutine chol_mo_from_chol_ao(cd_mo,cd_ao,n_mo,n_ao,n_cd,n_k,n_unique_k)
|
||||
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 MO integrals from 3-index AO integrals'
|
||||
print*,'providing 3-index CD MO integrals from 3-index CD AO integrals'
|
||||
|
||||
cd_mo = 0.d0
|
||||
|
||||
@ -315,7 +317,7 @@ subroutine chol_mo_from_chol_ao(cd_mo,cd_ao,n_mo,n_ao,n_cd,n_k,n_unique_k)
|
||||
)
|
||||
call wall_time(wall_2)
|
||||
call cpu_time(cpu_2)
|
||||
print*,' 3-idx MO provided'
|
||||
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),')'
|
||||
|
||||
|
||||
29
src/utils_complex/dump_cd_ints.irp.f
Normal file
29
src/utils_complex/dump_cd_ints.irp.f
Normal file
@ -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
|
||||
47
src/utils_complex/dump_cd_ksym.irp.f
Normal file
47
src/utils_complex/dump_cd_ksym.irp.f
Normal file
@ -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
|
||||
234
src/utils_complex/test_cd_ksym.irp.f
Normal file
234
src/utils_complex/test_cd_ksym.irp.f
Normal file
@ -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 <ij|kl> 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
|
||||
Loading…
Reference in New Issue
Block a user