LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-18 11:23:38 +01:00
Code Releases Activity

converter cleanup

Browse Source
This commit is contained in:
Kevin Gasperich 2020-03-11 13:48:35 -05:00
parent 8411167e90
commit f07bdee9cd
1 changed files with 225 additions and 148 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

373
src/utils_complex/MolPyscfToQPkpts.py
View File

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