mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-25 13:03:28 +01:00
Fixing converter for qmcpack and qp2 with pbc and excited states
This commit is contained in:
parent
4f9c398ffc
commit
671dd70ba8
@ -287,9 +287,12 @@ def convert_kpts(filename,qph5path,qmcpack=True):
|
|||||||
|
|
||||||
ezfio.set_qmcpack_qmc_pbc(qph5['qmcpack'].attrs['PBC'])
|
ezfio.set_qmcpack_qmc_pbc(qph5['qmcpack'].attrs['PBC'])
|
||||||
ezfio.set_qmcpack_qmc_cart(qph5['qmcpack'].attrs['cart'])
|
ezfio.set_qmcpack_qmc_cart(qph5['qmcpack'].attrs['cart'])
|
||||||
|
ezfio.set_qmcpack_qmc_pseudo(qph5['qmcpack'].attrs['Pseudo'])
|
||||||
ezfio.set_qmcpack_supertwist(qph5['qmcpack/Super_Twist'][()].tolist())
|
ezfio.set_qmcpack_supertwist(qph5['qmcpack/Super_Twist'][()].tolist())
|
||||||
ezfio.set_qmcpack_latticevectors(qph5['qmcpack/LatticeVectors'][()].tolist())
|
ezfio.set_qmcpack_latticevectors(qph5['qmcpack/LatticeVectors'][()].tolist())
|
||||||
|
ezfio.set_qmcpack_qmc_phase(qph5['qmcpack/qmc_phase'][()].tolist())
|
||||||
ezfio.set_qmcpack_qmc_mo_energy(qph5['qmcpack/eigenval'][()].tolist())
|
ezfio.set_qmcpack_qmc_mo_energy(qph5['qmcpack/eigenval'][()].tolist())
|
||||||
|
|
||||||
except AttributeError as err:
|
except AttributeError as err:
|
||||||
print("################################################")
|
print("################################################")
|
||||||
print("# ERROR: problem copying QMCPACK data to ezfio #")
|
print("# ERROR: problem copying QMCPACK data to ezfio #")
|
||||||
|
@ -398,7 +398,7 @@ def create_ezfio_stuff(dict_ezfio_cfg, config_or_default="config"):
|
|||||||
try:
|
try:
|
||||||
(begin, end) = list(map(str.strip, dim.split(":")))
|
(begin, end) = list(map(str.strip, dim.split(":")))
|
||||||
except ValueError:
|
except ValueError:
|
||||||
a_size_raw.append(dim)
|
a_size_raw.append(dim.strip())
|
||||||
else:
|
else:
|
||||||
if begin[0] == '-':
|
if begin[0] == '-':
|
||||||
a_size_raw.append("{0}+{1}+1".format(end, begin[1:]))
|
a_size_raw.append("{0}+{1}+1".format(end, begin[1:]))
|
||||||
|
@ -124,8 +124,7 @@ def makesq2(vlist,n1,n2):
|
|||||||
out[i] = tmp2.copy()
|
out[i] = tmp2.copy()
|
||||||
return out
|
return out
|
||||||
|
|
||||||
|
def get_supcellPhase(cell, kpts=[], kmesh=[]):
|
||||||
def get_phase(cell, kpts, kmesh=None):
|
|
||||||
'''
|
'''
|
||||||
The unitary transformation that transforms the supercell basis k-mesh
|
The unitary transformation that transforms the supercell basis k-mesh
|
||||||
adapted basis.
|
adapted basis.
|
||||||
@ -134,7 +133,7 @@ def get_phase(cell, kpts, kmesh=None):
|
|||||||
from pyscf import lib
|
from pyscf import lib
|
||||||
|
|
||||||
latt_vec = cell.lattice_vectors()
|
latt_vec = cell.lattice_vectors()
|
||||||
if kmesh is None:
|
if len(kmesh)== 0:
|
||||||
# Guess kmesh
|
# Guess kmesh
|
||||||
scaled_k = cell.get_scaled_kpts(kpts).round(8)
|
scaled_k = cell.get_scaled_kpts(kpts).round(8)
|
||||||
kmesh = (len(np.unique(scaled_k[:,0])),
|
kmesh = (len(np.unique(scaled_k[:,0])),
|
||||||
@ -153,58 +152,7 @@ def get_phase(cell, kpts, kmesh=None):
|
|||||||
|
|
||||||
# R_rel_mesh has to be construct exactly same to the Ts in super_cell function
|
# R_rel_mesh has to be construct exactly same to the Ts in super_cell function
|
||||||
scell = tools.super_cell(cell, kmesh)
|
scell = tools.super_cell(cell, kmesh)
|
||||||
return scell, phase
|
return scell,phase
|
||||||
|
|
||||||
def mo_k2gamma(cell, mo_energy, mo_coeff, kpts, kmesh=None):
|
|
||||||
'''
|
|
||||||
Transform MOs in Kpoints to the equivalents supercell
|
|
||||||
'''
|
|
||||||
from pyscf import lib
|
|
||||||
import scipy.linalg as la
|
|
||||||
scell, phase = get_phase(cell, kpts, kmesh)
|
|
||||||
|
|
||||||
E_g = np.hstack(mo_energy)
|
|
||||||
C_k = np.asarray(mo_coeff)
|
|
||||||
Nk, Nao, Nmo = C_k.shape
|
|
||||||
NR = phase.shape[0]
|
|
||||||
|
|
||||||
# Transform AO indices
|
|
||||||
C_gamma = np.einsum('Rk, kum -> Rukm', phase, C_k)
|
|
||||||
C_gamma = C_gamma.reshape(Nao*NR, Nk*Nmo)
|
|
||||||
|
|
||||||
E_sort_idx = np.argsort(E_g)
|
|
||||||
E_g = E_g[E_sort_idx]
|
|
||||||
C_gamma = C_gamma[:,E_sort_idx]
|
|
||||||
s = scell.pbc_intor('int1e_ovlp')
|
|
||||||
assert(abs(reduce(np.dot, (C_gamma.conj().T, s, C_gamma))
|
|
||||||
- np.eye(Nmo*Nk)).max() < 1e-7)
|
|
||||||
|
|
||||||
# Transform MO indices
|
|
||||||
E_k_degen = abs(E_g[1:] - E_g[:-1]).max() < 1e-5
|
|
||||||
if np.any(E_k_degen):
|
|
||||||
degen_mask = np.append(False, E_k_degen) | np.append(E_k_degen, False)
|
|
||||||
shift = min(E_g[degen_mask]) - .1
|
|
||||||
f = np.dot(C_gamma[:,degen_mask] * (E_g[degen_mask] - shift),
|
|
||||||
C_gamma[:,degen_mask].conj().T)
|
|
||||||
assert(abs(f.imag).max() < 1e-5)
|
|
||||||
|
|
||||||
e, na_orb = la.eigh(f.real, s, type=2)
|
|
||||||
C_gamma[:,degen_mask] = na_orb[:, e>0]
|
|
||||||
|
|
||||||
if abs(C_gamma.imag).max() < 1e-7:
|
|
||||||
print('!Warning Some complexe pollutions in MOs are present')
|
|
||||||
|
|
||||||
C_gamma = C_gamma.real
|
|
||||||
if abs(reduce(np.dot, (C_gamma.conj().T, s, C_gamma)) - np.eye(Nmo*Nk)).max() < 1e-7:
|
|
||||||
print('!Warning Some complexe pollutions in MOs are present')
|
|
||||||
|
|
||||||
s_k = cell.pbc_intor('int1e_ovlp', kpts=kpts)
|
|
||||||
# overlap between k-point unitcell and gamma-point supercell
|
|
||||||
s_k_g = np.einsum('kuv,Rk->kuRv', s_k, phase.conj()).reshape(Nk,Nao,NR*Nao)
|
|
||||||
# The unitary transformation from k-adapted orbitals to gamma-point orbitals
|
|
||||||
mo_phase = lib.einsum('kum,kuv,vi->kmi', C_k.conj(), s_k_g, C_gamma)
|
|
||||||
|
|
||||||
return mo_phase
|
|
||||||
|
|
||||||
def qp2rename():
|
def qp2rename():
|
||||||
import shutil
|
import shutil
|
||||||
@ -594,6 +542,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
|||||||
ne_threshold = int_threshold
|
ne_threshold = int_threshold
|
||||||
bielec_int_threshold = int_threshold
|
bielec_int_threshold = int_threshold
|
||||||
thresh_mono = int_threshold
|
thresh_mono = int_threshold
|
||||||
|
loc_cell,phase=get_supcellPhase(cell=cell,kmesh=kmesh,kpts=kpts)
|
||||||
|
|
||||||
|
|
||||||
# qph5path = 'qpdat.h5'
|
# qph5path = 'qpdat.h5'
|
||||||
@ -632,23 +581,39 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
|||||||
##########################################
|
##########################################
|
||||||
|
|
||||||
natom = cell.natm
|
natom = cell.natm
|
||||||
|
|
||||||
print('n_atom per kpt', natom)
|
print('n_atom per kpt', natom)
|
||||||
|
|
||||||
atom_xyz = mf.cell.atom_coords()
|
atom_xyz = mf.cell.atom_coords()
|
||||||
|
|
||||||
|
unit_bohr=1.0
|
||||||
if not(mf.cell.unit.startswith(('B','b','au','AU'))):
|
if not(mf.cell.unit.startswith(('B','b','au','AU'))):
|
||||||
from pyscf.data.nist import BOHR
|
from pyscf.data.nist import BOHR
|
||||||
atom_xyz /= BOHR # always convert to au
|
atom_xyz /= BOHR # always convert to au
|
||||||
|
unit_bohr=BOHR
|
||||||
|
|
||||||
with h5py.File(qph5path,'a') as qph5:
|
with h5py.File(qph5path,'a') as qph5:
|
||||||
qph5['nuclei'].attrs['kpt_num']=Nk
|
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)
|
if len(kpts)!=0:
|
||||||
atom_dset=qph5.create_dataset('nuclei/nucl_label',(natom,),dtype=strtype)
|
nbatom=loc_cell.natm
|
||||||
for i in range(natom):
|
qph5['nuclei'].attrs['nucl_num']=nbatom
|
||||||
atom_dset[i] = mf.cell.atom_pure_symbol(i)
|
MyPos=qph5.create_dataset('nuclei/nucl_coord',(nbatom,3),dtype="f8")
|
||||||
|
for x in range(nbatom):
|
||||||
|
MyPos[x:]=loc_cell.atom_coord(x)/unit_bohr
|
||||||
|
qph5.create_dataset('nuclei/nucl_charge',data=loc_cell.atom_charges())
|
||||||
|
strtype=h5py.special_dtype(vlen=str)
|
||||||
|
atom_dset=qph5.create_dataset('nuclei/nucl_label',(nbatom,),dtype=strtype)
|
||||||
|
for i in range(nbatom):
|
||||||
|
atom_dset[i] = loc_cell.atom_pure_symbol(i)
|
||||||
|
else:
|
||||||
|
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)
|
||||||
|
|
||||||
##########################################
|
##########################################
|
||||||
# #
|
# #
|
||||||
@ -725,6 +690,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
|||||||
qph5.create_dataset('qmcpack/qmc_lbas',data=qp_lbas)
|
qph5.create_dataset('qmcpack/qmc_lbas',data=qp_lbas)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
# with h5py.File(qph5path,'a') as qph5:
|
# with h5py.File(qph5path,'a') as qph5:
|
||||||
# qph5['mo_basis'].attrs['mo_num']=Nk*nmo
|
# qph5['mo_basis'].attrs['mo_num']=Nk*nmo
|
||||||
# qph5['ao_basis'].attrs['ao_num']=Nk*nao
|
# qph5['ao_basis'].attrs['ao_num']=Nk*nao
|
||||||
@ -890,90 +856,17 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
|
|||||||
if len(kpts)== 0:
|
if len(kpts)== 0:
|
||||||
sp_twist=[0.0,0.0,0.0]
|
sp_twist=[0.0,0.0,0.0]
|
||||||
|
|
||||||
|
|
||||||
with h5py.File(qph5path,'a') as qph5:
|
with h5py.File(qph5path,'a') as qph5:
|
||||||
qph5['qmcpack'].attrs['PBC']=True
|
qph5['qmcpack'].attrs['PBC']=True
|
||||||
qph5['qmcpack'].attrs['cart']=cell.cart
|
qph5['qmcpack'].attrs['cart']=cell.cart
|
||||||
|
qph5['qmcpack'].attrs['Pseudo']=cell.has_ecp()
|
||||||
qph5.create_dataset('qmcpack/Super_Twist',(1,3),dtype="f8",data=sp_twist)
|
qph5.create_dataset('qmcpack/Super_Twist',(1,3),dtype="f8",data=sp_twist)
|
||||||
qph5.create_dataset('qmcpack/LatticeVectors',(3,3),dtype="f8",data=cell.lattice_vectors().T)
|
if len(kpts)!= 0:
|
||||||
|
qph5.create_dataset('qmcpack/LatticeVectors',(3,3),dtype="f8",data=loc_cell.lattice_vectors())
|
||||||
|
else:
|
||||||
|
qph5.create_dataset('qmcpack/LatticeVectors',(3,3),dtype="f8",data=cell.lattice_vectors())
|
||||||
qph5.create_dataset('qmcpack/eigenval',(1,Nk*nmo),dtype="f8",data=mf.mo_energy)
|
qph5.create_dataset('qmcpack/eigenval',(1,Nk*nmo),dtype="f8",data=mf.mo_energy)
|
||||||
|
qph5.create_dataset('qmcpack/qmc_phase',data=phase.view(dtype=float))
|
||||||
|
|
||||||
# ##########################################
|
|
||||||
# # #
|
|
||||||
# # ECP #
|
|
||||||
# # #
|
|
||||||
# ##########################################
|
|
||||||
#
|
|
||||||
# if (cell.has_ecp()):
|
|
||||||
# #atsymb = [mol.atom_pure_symbol(i) for i in range(natom)]
|
|
||||||
# #pyecp = mol._ecp
|
|
||||||
# ## nelec to remove for each atom
|
|
||||||
# #nuc_z_remov = [pyecp[i][0] for i in atsymb]
|
|
||||||
# #nl_per_atom = [len(pyecp[i][1]) for i in atsymb]
|
|
||||||
# ## list of l-values for channels of each atom
|
|
||||||
# #ecp_l = [[pyecp[i][1][j][0] for j in range(len(pyecp[i][1]))] for i in atsymb]
|
|
||||||
# ## list of [exp,coef] for each channel (r**0,1,2,3,4,5,)
|
|
||||||
# #ecp_ac = [[pyecp[i][1][j][1] for j in range(len(pyecp[i][1]))] for i in atsymb]
|
|
||||||
# pyecp = [cell._ecp[cell.atom_pure_symbol(i)] for i in range(natom)]
|
|
||||||
# nzrmv=[0]*natom
|
|
||||||
# lmax=0
|
|
||||||
# klocmax=0
|
|
||||||
# knlmax=0
|
|
||||||
# for i,(nz,dat) in enumerate(pyecp):
|
|
||||||
# nzrmv[i]=nz
|
|
||||||
# for lval,ac in dat:
|
|
||||||
# if (lval==-1):
|
|
||||||
# klocmax=max(sum(len(j) for j in ac),klocmax)
|
|
||||||
# else:
|
|
||||||
# lmax=max(lval,lmax)
|
|
||||||
# knlmax=max(sum(len(j) for j in ac),knlmax)
|
|
||||||
# #psd_nk = np.zeros((natom,klocmax),dtype=int)
|
|
||||||
# #psd_vk = np.zeros((natom,klocmax),dtype=float)
|
|
||||||
# #psd_dzk = np.zeros((natom,klocmax),dtype=float)
|
|
||||||
# #psd_nkl = np.zeros((natom,knlmax,lmax+1),dtype=int)
|
|
||||||
# #psd_vkl = np.zeros((natom,knlmax,lmax+1),dtype=float)
|
|
||||||
# #psd_dzkl = np.zeros((natom,knlmax,lmax+1),dtype=float)
|
|
||||||
# klnlmax=max(klocmax,knlmax)
|
|
||||||
# psd_n = np.zeros((lmax+2,klnlmax,natom),dtype=int)
|
|
||||||
# psd_v = np.zeros((lmax+2,klnlmax,natom),dtype=float)
|
|
||||||
# psd_dz = np.zeros((lmax+2,klnlmax,natom),dtype=float)
|
|
||||||
# for i,(_,dat) in enumerate(pyecp):
|
|
||||||
# for lval,ac in dat:
|
|
||||||
# count=0
|
|
||||||
# for ri,aici in enumerate(ac):
|
|
||||||
# for ai,ci in aici:
|
|
||||||
# psd_n[lval+1,count,i] = ri-2
|
|
||||||
# psd_v[lval+1,count,i] = ci
|
|
||||||
# psd_dz[lval+1,count,i] = ai
|
|
||||||
# count += 1
|
|
||||||
# psd_nk = psd_n[0,:klocmax]
|
|
||||||
# psd_vk = psd_v[0,:klocmax]
|
|
||||||
# psd_dzk = psd_dz[0,:klocmax]
|
|
||||||
# psd_nkl = psd_n[1:,:knlmax]
|
|
||||||
# psd_vkl = psd_v[1:,:knlmax]
|
|
||||||
# psd_dzkl = psd_dz[1:,:knlmax]
|
|
||||||
# with h5py.File(qph5path,'a') as qph5:
|
|
||||||
# qph5['pseudo'].attrs['do_pseudo']=True
|
|
||||||
# qph5['pseudo'].attrs['pseudo_lmax']=lmax
|
|
||||||
# qph5['pseudo'].attrs['pseudo_klocmax']=klocmax
|
|
||||||
# qph5['pseudo'].attrs['pseudo_kmax']=knlmax
|
|
||||||
# qph5.create_dataset('pseudo/nucl_charge_remove',data=nzrmv)
|
|
||||||
# qph5.create_dataset('pseudo/pseudo_n_k',data=psd_nk)
|
|
||||||
# qph5.create_dataset('pseudo/pseudo_n_kl',data=psd_nkl)
|
|
||||||
# qph5.create_dataset('pseudo/pseudo_v_k',data=psd_vk)
|
|
||||||
# qph5.create_dataset('pseudo/pseudo_v_kl',data=psd_vkl)
|
|
||||||
# qph5.create_dataset('pseudo/pseudo_dz_k',data=psd_dzk)
|
|
||||||
# qph5.create_dataset('pseudo/pseudo_dz_kl',data=psd_dzkl)
|
|
||||||
#
|
|
||||||
# ## nelec to remove for each atom
|
|
||||||
# #nuc_z_remov = [i[0] for i in pyecp]
|
|
||||||
# #nl_per_atom = [len(i[1]) for i in pyecp]
|
|
||||||
# ## list of l-values for channels of each atom
|
|
||||||
# #ecp_l = [[ j[0] for j in i[1] ] for i in pyecp]
|
|
||||||
# #lmax = max(map(max,ecp_l))
|
|
||||||
# ## list of [exp,coef] for each channel (r**0,1,2,3,4,5,)
|
|
||||||
# #ecp_ac = [[ j[1] for j in i[1] ] for i in pyecp]
|
|
||||||
|
|
||||||
return
|
return
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user