LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-27 18:54:42 +02:00
Code Releases Activity

Merge pull request #19 from kgasperich/cd-pbc-conv-patch

Browse Source
This commit is contained in:
Kevin Gasperich 2022-10-06 12:25:16 -05:00 committed by GitHub
commit 54ce1c24f5
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
17 changed files with 2210 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -3,7 +3,6 @@ quantum_package_static.tar.gz
build.ninja build.ninja
.ninja_log .ninja_log
.ninja_deps .ninja_deps
bin/
lib/ lib/
config/qp_create_ninja.pickle config/qp_create_ninja.pickle
src/*/.gitignore src/*/.gitignore

720
bin/qp_convert_h5_to_ezfio
View File

@ -3,12 +3,13 @@
convert hdf5 output (e.g. from PySCF) to ezfio convert hdf5 output (e.g. from PySCF) to ezfio
Usage: Usage:
qp_convert_h5_to_ezfio [--noqmc] [-o EZFIO_DIR] FILE qp_convert_h5_to_ezfio [--noqmc] [--rmg] [-o EZFIO_DIR] FILE
Options: Options:
-o --output=EZFIO_DIR Produced directory -o --output=EZFIO_DIR Produced directory
by default is FILE.ezfio by default is FILE.ezfio
--noqmc don't include basis, cell, etc. for QMCPACK --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 from ezfio import ezfio
@ -20,6 +21,104 @@ from docopt import docopt
import gzip import gzip
#fname = sys.argv[1] #fname = sys.argv[1]
#qph5name = sys.argv[2] #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)
kcon_p = np.zeros((nk,nk,nk),dtype=int)
for i in range(nk):
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
def get_full_path(file_path): def get_full_path(file_path):
file_path = os.path.expanduser(file_path) file_path = os.path.expanduser(file_path)
@ -27,6 +126,19 @@ def get_full_path(file_path):
# file_path = os.path.abspath(file_path) # file_path = os.path.abspath(file_path)
return file_path return file_path
def make_reim_identity_kblocks(nk,nm,na=None):
if na is None:
na = nm
single_block = np.eye(nm, na, dtype=np.complex128).view(dtype=np.float64).reshape((nm, na, 2))
kblocks = np.tile(single_block,[nk, 1, 1, 1])
return kblocks
def make_reim_identity_block_diag(nk,nm,na=None):
from scipy.linalg import block_diag
kblocks = make_reim_identity_kblocks(nk,nm,na).view(dtype=np.complex128).squeeze()
kblockdiag = block_diag(*kblocks).view(dtype=np.float64).reshape((nk*nm,nk*na,2))
print(f'kblockdiag.shape = {kblockdiag.shape}')
return kblockdiag
def flatten(l): def flatten(l):
res = [] res = []
@ -173,6 +285,592 @@ def convert_mol(filename,qph5path):
return return
def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
import json
from scipy.linalg import block_diag
dump_fci, dump_cd, dump_fci2 = (False, False, False)
ezfio.set_file(filename)
ezfio.set_nuclei_is_complex(True)
# Dummy atom since AFQMC h5 has no atom information
#nucl_num = 1
#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:
kpt_num = qph5['Hamiltonian/KPoints'][()].shape[0]
ham_dims = qph5['Hamiltonian/dims'][()]
NMOPerKP = qph5['Hamiltonian/NMOPerKP'][()]
_, _, kpt_num, orb_num, elec_alpha_num_tot, elec_beta_num_tot, _, nchol_maybe = ham_dims
#for now, all kpts must have same number of MOs
for nmoi in NMOPerKP:
if nmoi != NMOPerKP[0]:
print("ERROR: all KPs must have same number of MOs")
raise ValueError
#TODO: fix na, nb in rmg
assert(elec_alpha_num_tot % kpt_num == 0)
assert(elec_beta_num_tot % kpt_num == 0)
elec_alpha_num_per_kpt = elec_alpha_num_tot // kpt_num
elec_beta_num_per_kpt = elec_beta_num_tot // kpt_num
#elec_alpha_num_per_kpt = qph5['Hamiltonian']['dims'][4]
#elec_beta_num_per_kpt = qph5['Hamiltonian']['dims'][5]
#orb_num = qph5['Hamiltonian']['dims'][3]
#try:
# is_ao = json.loads(qph5['metadata'][()].decode("utf-8").replace("'",'"'))['ortho_ao']
# if is_ao:
# ao_num = orb_num
# elif is_ao ==False:
# mo_num = orb_num
# else:
# raise ValueError('Problem with ortho_ao key in metadata')
#except:
# raise UnicodeDecodeError('metadata not correctly parsed from HDF5 file')
ezfio.set_nuclei_kpt_num(kpt_num)
kpt_pair_num = (kpt_num*kpt_num + kpt_num)//2
ezfio.set_nuclei_kpt_pair_num(kpt_pair_num)
# don't multiply nuclei by kpt_num
# work in k-space, not in equivalent supercell
#nucl_num_per_kpt = nucl_num
#ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)
# these are totals (kpt_num * num_per_kpt)
# need to change if we want to truncate orbital space within pyscf
#if is_ao:
# ao_num = orb_num*kpt_num
#TODO: fix this?
ao_num_tot = orb_num
ao_num_per_kpt = NMOPerKP[0]
mo_num_tot = orb_num
mo_num_per_kpt = NMOPerKP[0]
#mo_num_per_kpt = ao_num_per_kpt
ezfio.set_ao_basis_ao_num(ao_num_per_kpt * kpt_num)
ezfio.set_mo_basis_mo_num(mo_num_per_kpt * kpt_num)
ezfio.set_ao_basis_ao_num_per_kpt(ao_num_per_kpt)
ezfio.set_mo_basis_mo_num_per_kpt(mo_num_per_kpt)
ezfio.electrons_elec_alpha_num = elec_alpha_num_per_kpt * kpt_num
ezfio.electrons_elec_beta_num = elec_beta_num_per_kpt * kpt_num
##########################################
# #
# Basis #
# #
##########################################
#TODO
nucl_num = 1
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 )
nucl_num_per_kpt = nucl_num
ezfio.set_nuclei_nucl_num(nucl_num_per_kpt)
ezfio.set_nuclei_io_kpt_symm('Read')
ezfio.set_ao_basis_ao_basis("dummy basis")
#nucleus on which each AO is centered
ao_nucl = [1 for i in range(ao_num_per_kpt)]*kpt_num
ezfio.set_ao_basis_ao_nucl(ao_nucl)
#Just need one (can clean this up later)
ao_prim_num_max = 5
d = [ [0] *ao_prim_num_max]*ao_num_tot
ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num_tot)
ezfio.set_ao_basis_ao_power(d)
ezfio.set_ao_basis_ao_coef(d)
ezfio.set_ao_basis_ao_expo(d)
##########################################
# #
# MOCoeff #
# #
##########################################
#TODO
#coef_per_kpt = np.eye(mo_num_per_kpt, ao_num_per_kpt, dtype=np.complex128).view(dtype=np.float64).reshape((mo_num_per_kpt, ao_num_per_kpt, 2))
#mo_coef_kpts = np.tile(coef_per_kpt,[kpt_num, 1, 1, 1])
#qph5.create_dataset('mo_basis/mo_coef_kpts',data=make_reim_identity_kblocks(kpt_num, mo_num_per_kpt, ao_num_per_kpt))
ezfio.set_mo_basis_mo_coef_kpts(make_reim_identity_kblocks(kpt_num, mo_num_per_kpt, ao_num_per_kpt))
ezfio.set_mo_basis_mo_coef_complex(make_reim_identity_block_diag(kpt_num, mo_num_per_kpt, ao_num_per_kpt))
##########################################
# #
# Integrals Mono #
# #
##########################################
with h5py.File(qph5path,'r') as qph5:
# we don't have separate kinetic, nuc-elec, pseudo 1e ints, so just combine in nuc-elec and set rest to zero
mono_ints_tot = np.zeros((kpt_num,ao_num_per_kpt,ao_num_per_kpt,2),dtype=np.float64)
for i in range(kpt_num):
mono_ints_tot[i] = qph5[f'Hamiltonian/H1_kp{i}'][()]
ovlp_ao_reim = make_reim_identity_kblocks(kpt_num,ao_num_per_kpt,ao_num_per_kpt)
kin_ao_reim = np.zeros((kpt_num,ao_num_per_kpt,ao_num_per_kpt,2),dtype=np.float64)
ne_ao_reim = mono_ints_tot
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 is_ao:
kin_ao_reim=
ovlp_ao_reim=
ne_ao_reim=
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')
else:
kin_mo_reim=
ovlp_mo_reim=
ne_mo_reim=
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_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')
"""
##########################################
# #
# k-points #
# #
##########################################
#TODO
with h5py.File(qph5path,'r') as qph5:
#kconserv = qph5['nuclei/kconserv'][()].tolist()
minusk = qph5['Hamiltonian']['MinusK'][:]+1
QKTok2 = qph5['Hamiltonian']['QKTok2'][:]+1
#TODO: change this after rmg is fixed
#minusk = QKTok2[:,0]
kconserv = kconserv_p_from_qkk2_mk(QKTok2-1,minusk-1)+1
unique_kpt_num = len(qph5['Hamiltonian']['KPFactorized'])
unique_k_idx = []
for i in qph5['Hamiltonian']['KPFactorized'].keys():
unique_k_idx.append(int(i[1:])+1)
unique_k_idx.sort()
kpt_sparse_map = np.zeros(kpt_num,dtype=int)
isparse=0
#TODO: make robust: this assumes that for each pair, the one with data has a lower index
for i in range(kpt_num):
if i+1 in unique_k_idx:
kpt_sparse_map[i] = isparse+1
isparse += 1
else:
kpt_sparse_map[i] = -kpt_sparse_map[minusk[i]-1]
ezfio.set_nuclei_kconserv(kconserv.transpose(2,1,0))
ezfio.set_nuclei_io_kconserv('Read')
ezfio.set_nuclei_minusk(minusk)
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
# unique_kpt_num
# io_kpt_symm
##########################################
# #
# Integrals Bi #
# #
##########################################
# should this be in ao_basis? ao_two_e_ints?
with h5py.File(qph5path,'r') as qph5:
nchol_per_kpt_all = qph5['Hamiltonian']['NCholPerKP'][:]
print(f'nchol_per_kpt_full = {nchol_per_kpt_all}')
#nchol_per_kpt = nchol_per_kpt_all[nchol_per_kpt_all != 0]
nchol_per_kpt = nchol_per_kpt_all[np.array(unique_k_idx,dtype=int)-1]
print(f'nchol_per_kpt = {nchol_per_kpt}')
print(f'unique_k_idx = {unique_k_idx}')
#for i in range(kpt_num):
# if i+1 in unique_k_idx:
# print('* ',i,nchol_per_kpt_all[i])
# else:
# print(' ',i,nchol_per_kpt_all[i])
nchol_per_kpt_max = max(nchol_per_kpt)
ezfio.set_ao_two_e_ints_chol_num(nchol_per_kpt)
ezfio.set_ao_two_e_ints_chol_num_max(nchol_per_kpt_max)
if is_ao:
#ao_num_per_kpt = ao_num//kpt_num
ezfio.set_ao_two_e_ints_io_chol_ao_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 = []
L_all = np.zeros((unique_kpt_num, kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt_max,2),dtype=np.float64)
print(f'kpt_sparse_map = {kpt_sparse_map}')
print(f'unique_k_idx-1 = {np.array(unique_k_idx)-1}')
for i in range(unique_kpt_num):
ki = unique_k_idx[i]-1
#print(i, ki)
L_i = qph5[f'Hamiltonian/KPFactorized/L{ki}'][()].reshape((kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[kpt_sparse_map[ki]-1], 2))
#L.reshape(kpt_num, ao_num_per_kpt, ao_num_per_kpt, nchol_per_kpt[i], 2)
#L = np.einsum("ijklm->ilkjm", 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): def convert_kpts(filename,qph5path,qmcpack=True):
ezfio.set_file(filename) ezfio.set_file(filename)
ezfio.set_nuclei_is_complex(True) ezfio.set_nuclei_is_complex(True)
@ -634,6 +1332,7 @@ def convert_cplx(filename,qph5path):
ezfio.set_nuclei_kconserv(kconserv) ezfio.set_nuclei_kconserv(kconserv)
ezfio.set_nuclei_io_kconserv('Read') ezfio.set_nuclei_io_kconserv('Read')
########################################## ##########################################
# # # #
# Integrals Bi # # Integrals Bi #
@ -671,8 +1370,15 @@ def convert_cplx(filename,qph5path):
if __name__ == '__main__': if __name__ == '__main__':
ARGUMENTS = docopt(__doc__) 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']) FILE = get_full_path(ARGUMENTS['FILE'])
qmcpack = True qmcpack = True
rmg = False
if ARGUMENTS["--output"]: if ARGUMENTS["--output"]:
EZFIO_FILE = get_full_path(ARGUMENTS["--output"]) EZFIO_FILE = get_full_path(ARGUMENTS["--output"])
else: else:
@ -680,10 +1386,19 @@ if __name__ == '__main__':
if ARGUMENTS["--noqmc"]: if ARGUMENTS["--noqmc"]:
qmcpack = False qmcpack = False
if ARGUMENTS["--rmg"]:
rmg = True
with h5py.File(FILE,'r') as qph5: with h5py.File(FILE,'r') as qph5:
try:
do_kpts = ('kconserv' in qph5['nuclei'].keys()) do_kpts = ('kconserv' in qph5['nuclei'].keys())
if (do_kpts): except:
do_kpts = False
if (do_kpts or rmg):
print("converting HDF5 to EZFIO for periodic system") print("converting HDF5 to EZFIO for periodic system")
if rmg:
print("Using RMG and AFQMC h5")
convert_kpts_cd(EZFIO_FILE,FILE,qmcpack)
else:
convert_kpts(EZFIO_FILE,FILE,qmcpack) convert_kpts(EZFIO_FILE,FILE,qmcpack)
else: else:
print("converting HDF5 to EZFIO for molecular system") print("converting HDF5 to EZFIO for molecular system")
@ -697,4 +1412,3 @@ if __name__ == '__main__':
# #
#to be sure your MOs will be orthogonal, which is not the case when #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).""") #the MOs are read from output files (not enough precision in output).""")

23
src/ao_two_e_ints/EZFIO.cfg
View File

@ -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) 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 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

261
src/ao_two_e_ints/cd_ao_ints.irp.f Normal file
View File

@ -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

5
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -358,6 +358,11 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ]
print*, 'AO integrals provided from 3-index ao ints (periodic)' print*, 'AO integrals provided from 3-index ao ints (periodic)'
ao_two_e_integrals_in_map = .True. ao_two_e_integrals_in_map = .True.
return 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 else
print*,'calculation of periodic AOs not implemented' print*,'calculation of periodic AOs not implemented'
stop -1 stop -1

19
src/hartree_fock/hf_energy.irp.f
View File

@ -13,19 +13,22 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, hf_energy] BEGIN_PROVIDER [ double precision, hf_energy]
&BEGIN_PROVIDER [ double precision, hf_two_electron_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] &BEGIN_PROVIDER [ double precision, hf_one_electron_energy]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components.
END_DOC END_DOC
integer :: i,j,k integer :: i,j,k,jk
hf_energy = nuclear_repulsion hf_energy = nuclear_repulsion
hf_two_electron_energy = 0.d0 hf_two_electron_energy = 0.d0
hf_two_electron_energy_jk = 0.d0
hf_one_electron_energy = 0.d0 hf_one_electron_energy = 0.d0
if (is_complex) then 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_1e_tmp = (0.d0,0.d0)
hf_2e_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 k=1,kpt_num
do j=1,ao_num_per_kpt do j=1,ao_num_per_kpt
do i=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) ) +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) & 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) ) + 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
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 if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then
print*,'HF_2e energy should be real:',irp_here print*,'HF_2e energy should be real:',irp_here
stop -1 stop -1

3
src/hartree_fock/print_e_scf.irp.f
View File

@ -15,6 +15,9 @@ subroutine run
print*,hf_one_electron_energy print*,hf_one_electron_energy
print*,hf_two_electron_energy print*,hf_two_electron_energy
print*,hf_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 end

12
src/mo_two_e_ints/EZFIO.cfg
View File

@ -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) 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 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

325
src/mo_two_e_ints/cd_mo_ints.irp.f Normal file
View File

@ -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

6
src/mo_two_e_ints/mo_bi_integrals.irp.f
View File

@ -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_single
call mo_map_fill_from_df_dot call mo_map_fill_from_df_dot
return 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 else
PROVIDE ao_two_e_integrals_in_map PROVIDE ao_two_e_integrals_in_map
endif endif

29
src/nuclei/EZFIO.cfg
View File

@ -60,3 +60,32 @@ type: integer
doc: array containing information about k-point symmetry doc: array containing information about k-point symmetry
size: (nuclei.kpt_num,nuclei.kpt_num,nuclei.kpt_num) size: (nuclei.kpt_num,nuclei.kpt_num,nuclei.kpt_num)
interface: ezfio 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

101
src/nuclei/kconserv_cplx.irp.f
View File

@ -21,8 +21,9 @@ BEGIN_PROVIDER [integer, kconserv, (kpt_num,kpt_num,kpt_num)]
call ezfio_get_nuclei_kconserv(kconserv) call ezfio_get_nuclei_kconserv(kconserv)
print *, 'kconserv read from disk' print *, 'kconserv read from disk'
else else
print*,'kconserv must be provided' call set_kconserv(kconserv)
stop -1 !print*,'kconserv must be provided'
!stop -1
endif endif
if (write_kconserv) then if (write_kconserv) then
call ezfio_set_nuclei_kconserv(kconserv) call ezfio_set_nuclei_kconserv(kconserv)
@ -30,6 +31,86 @@ BEGIN_PROVIDER [integer, kconserv, (kpt_num,kpt_num,kpt_num)]
endif endif
END_PROVIDER 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) subroutine double_allowed_kpts(kh1,kh2,kp1,kp2,is_allowed)
implicit none implicit none
integer, intent(in) :: kh1,kh2,kp1,kp2 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) is_allowed = (kconserv(kh1,kh2,kp1) == kp2)
end subroutine 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

259
src/scf_utils/fock_matrix_cplx.irp.f
View File

@ -540,6 +540,265 @@ BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_kpts, (ao_num_per_kpt, ao_num_per_kp
endif endif
enddo enddo
endif endif
END_PROVIDER
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 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_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
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, &
!$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)*(<ij|kl>)
!G_b(i,k) += D_{ab}(l,j)*(<ij|kl>)
!G_a(i,l) -= D_a (k,j)*(<ij|kl>)
!G_b(i,l) -= D_b (k,j)*(<ij|kl>)
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)*(<ij|kl>)
!G_b(i,k) += D_{ab}(l,j)*(<ij|kl>)
!G_a(i,l) -= D_a (k,j)*(<ij|kl>)
!G_b(i,l) -= D_b (k,j)*(<ij|kl>)
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 END_PROVIDER

24
src/utils_complex/dump_ao_2e_cplx.irp.f
View File

@ -15,19 +15,21 @@ subroutine run
do k=1,ao_num do k=1,ao_num
do l=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) 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
enddo enddo
enddo enddo
print*,'map1' !print*,'map1'
do i=0,ao_integrals_map%map_size !do i=0,ao_integrals_map%map_size
print*,i,ao_integrals_map%map(i)%value(:) ! print*,i,ao_integrals_map%map(i)%value(:)
print*,i,ao_integrals_map%map(i)%key(:) ! print*,i,ao_integrals_map%map(i)%key(:)
enddo !enddo
print*,'map2' !print*,'map2'
do i=0,ao_integrals_map_2%map_size !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)%value(:)
print*,i,ao_integrals_map_2%map(i)%key(:) ! print*,i,ao_integrals_map_2%map(i)%key(:)
enddo !enddo
end end

29
src/utils_complex/dump_cd_ints.irp.f Normal file
View 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
View 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
View 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