From 73cb1e9cd786d73c01b61715fbd427cf563354be Mon Sep 17 00:00:00 2001 From: Woochang Kim Date: Wed, 19 Jul 2023 21:30:41 -0700 Subject: [PATCH] update hartree --- WS2/LDA_SG15/1-scf/read_pw2bgw_grid.py | 49 ++++++ WS2/LDA_SG15/3.0-libxc/libxc_test.py | 29 ++-- WS2/LDA_SG15/4.0-hartree/compute_hartree.py | 36 ++-- WS2/LDA_SG15/4.0-hartree/density2Vh.py | 172 ++++++++++++++++++++ WS2/LDA_SG15/4.0-hartree/fft_mod.py | 65 +++++++- WS2/LDA_SG15/4.0-hartree/fft_test.py | 165 +++++++++++++++++++ WS2/LDA_SG15/4.0-hartree/fft_tutorial.py | 35 ++++ 7 files changed, 520 insertions(+), 31 deletions(-) create mode 100644 WS2/LDA_SG15/1-scf/read_pw2bgw_grid.py create mode 100644 WS2/LDA_SG15/4.0-hartree/density2Vh.py create mode 100644 WS2/LDA_SG15/4.0-hartree/fft_test.py create mode 100644 WS2/LDA_SG15/4.0-hartree/fft_tutorial.py diff --git a/WS2/LDA_SG15/1-scf/read_pw2bgw_grid.py b/WS2/LDA_SG15/1-scf/read_pw2bgw_grid.py new file mode 100644 index 0000000..e4a8a84 --- /dev/null +++ b/WS2/LDA_SG15/1-scf/read_pw2bgw_grid.py @@ -0,0 +1,49 @@ +import h5py +import numpy as np + +f = h5py.File('./WS2.save/charge-density.hdf5','r') +print(list(f.keys())) +print() +['MillerIndices', 'rhotot_g'] +igvec = f['MillerIndices'][:] +rhog = f['rhotot_g'][:] +ngm_g = f.attrs['ngm_g'] +print(ngm_g) +print(igvec.dtype) +print(igvec.shape) +print(rhog.dtype) +print(rhog.shape) +print(rhog[ngm_g]) + +print(rhog[0]*1449.7341) +print(rhog[1]*1449.7341) +rhog = rhog[::2] + 1j*rhog[1::2] +print(rhog.shape) +np.save('rhog', rhog) +np.save('igvec', igvec) +##data = np.loadtxt('./VXC.txt', skiprows=13, dtype=np.complex128) +##print(data) +#f = open('./pw2bgw.save/VXC.txt', 'r') +#f.readline() +##nsf, ng_g, ntran, cell_symmetry, nat, ecutrho = np.int32(f.readline().split()) +#f.readline() +##n1, n2, n3 = no.int32(f.readline().split()) +#f.readline() +#f.readline() +#f.readline() +#f.readline() +#f.readline() +#f.readline() +#f.readline() +#ng = np.int32(f.readline().split()[0]) +#gvec = np.int32(f.readline().split()).reshape((ng,3)) +#print(gvec[:5]) +#f.readline() +#ng = np.int32(f.readline().split()[0]) +#data = f.readline().split() +#f.close() +#data = np.array([np.float64(eval(str)) for str in data]).view(np.complex128).reshape(ng) +#print(data[0:5]) +#print(data.shape) +#np.save('vxc_gvec',gvec) +#np.save('vxc_pw2bgw',data) diff --git a/WS2/LDA_SG15/3.0-libxc/libxc_test.py b/WS2/LDA_SG15/3.0-libxc/libxc_test.py index a96574a..44ea4ae 100644 --- a/WS2/LDA_SG15/3.0-libxc/libxc_test.py +++ b/WS2/LDA_SG15/3.0-libxc/libxc_test.py @@ -3,27 +3,29 @@ from fft_mod import put_FFTbox2, ifft_g2r from pymatgen.io.cube import Cube from write_cube_mod import write_cube +bohr2ang = 0.529177210903 # Load data def vxc_libxc(): n1, n2, n3 = [18, 18, 135] #gvec_rho = np.loadtxt('../2.1-wfn/rho_gvec.txt',dtype=np.int32) #rho_pw2bgw = np.load('../2.1-wfn/rho_pw2bgw.npy') - #rho_FFTbox = put_FFTbox2(rho_pw2bgw, gvec_rho, [n1,n2,n3], noncolin=False) - #rho = ifft_g2r(rho_FFTbox) + gvec_rho = np.load('../1-scf/igvec.npy') + rho_pw2bgw = np.load('../1-scf/rhog.npy')*18*18*135 + rho_FFTbox = put_FFTbox2(rho_pw2bgw, gvec_rho, [n1,n2,n3], noncolin=False) + rho = ifft_g2r(rho_FFTbox) #rho_pp = Cube("../1.1-pp/Rho.cube") #rho = rho_pp.data - rho_pw2bgw = np.load("../2.1-wfn/rhor_pw2bgw.npy") - rho = rho_pw2bgw[0,:,:,:] + #rho_pw2bgw = np.load("../2.1-wfn/rhor_pw2bgw.npy") + #rho = rho_pw2bgw[0,:,:,:] # Create input - const = 1 #const = 214.82791025888773*18*18*135 #const = 18*18*135 inp = {} rho_flatten = rho.flatten() #print(rho - np.reshape(rho_flatten, (18,18, 135))) - inp["rho"] = const*rho_flatten + inp["rho"] = rho_flatten # Build functional func_c = pylibxc.LibXCFunctional("LDA_C_PZ", "unpolarized") @@ -47,14 +49,21 @@ def vxc_libxc(): #vxc = zk_x + zk_c + inp["rho"]*vrho_x + inp["rho"]*vrho_c vxc = vrho_x + vrho_c + zk = zk_x + zk_c vxc_FFTbox = np.reshape(vxc,(18,18,135)) + zk_FFTbox = np.reshape(zk,(18,18,135)) + Exc = np.sum(zk_FFTbox*rho) np.save('vxc_FFTbox', vxc_FFTbox) - return vxc_FFTbox + return vxc_FFTbox, Exc -vxc_gen = vxc_libxc() +vxc_gen, Exc = vxc_libxc() print(vxc_gen.shape) +n1, n2, n3 = [18, 18, 135] +nfft = n1*n2*n3 vxc_pp = Cube("../1.1-pp/Vxc.cube") +vol = vxc_pp.volume/(bohr2ang**3) # in Bohr^3 +vol = 1449.7341 vxc_pp = vxc_pp.data #vxc_pp = np.load('../2.1-wfn/vxcr_pw2bgw.npy') #vxc_pp = vxc_pp[0,:,:,:] @@ -67,11 +76,9 @@ def vxc_libxc(): #write_cube('ratio.cube',ratio,struct) #print(ratio) #print('\n') -print(np.max(ratio)) -print(np.min(ratio)) -print(np.average((ratio))) diff = np.abs(2*vxc_gen - vxc_pp) print(np.max(diff)) print(np.min(diff)) print(np.average(diff)) #print(np.sort(ratio.flat)[-1000:]) +print('Exc*2, ',2*Exc*(vol/nfft)) diff --git a/WS2/LDA_SG15/4.0-hartree/compute_hartree.py b/WS2/LDA_SG15/4.0-hartree/compute_hartree.py index 96932c7..c73e25e 100644 --- a/WS2/LDA_SG15/4.0-hartree/compute_hartree.py +++ b/WS2/LDA_SG15/4.0-hartree/compute_hartree.py @@ -1,39 +1,53 @@ import numpy as np -from fft_mod import put_FFTbox2, ifft_g2r +from fft_mod import put_FFTbox2, ifft_g2r, fft_r2g, flat_FFTbox from pymatgen.io.cube import Cube -bohr2ang = 0.529177 +bohr2ang = 0.529177249 def main(): - rho_g = np.load("../2.1-wfn/rho_pw2bgw.npy") - igvec = np.load("../2.1-wfn/rho_gvec.npy") + #rho_g = np.load("../2.1-wfn/rho_pw2bgw.npy") + #igvec = np.load("../2.1-wfn/rho_gvec.npy") + igvec = np.load('../1-scf/igvec.npy') + rho_g = np.load('../1-scf/rhog.npy') n1, n2, n3 = [18, 18, 135] nfft = n1*n2*n3 vbh = Cube("../1.1-pp/Vh.cube") #just for struct struct = vbh.structure vol = vbh.volume/(bohr2ang**3) # in Bohr^3 + vol = 1449.7341 cell = struct.lattice._matrix*(1/bohr2ang) # in Bohr bcell = 2*np.pi*np.linalg.inv(cell).T # in Bohr^-1 # debug - rho_g /= (vol/nfft) + #rho_g /= (vol/nfft) + rho_g *= (nfft) rho_g_FFTbox = put_FFTbox2(rho_g, igvec, [n1, n2, n3], noncolin=False) - rho_r = ifft_g2r(rho_g_FFTbox) - print('np.sum(rho_r)*(vol/nfft)',np.sum(rho_r)*(vol/nfft)) + rho_r_FFTbox = ifft_g2r(rho_g_FFTbox) + #rho_r_FFTbox = np.zeros((1,n1,n2,n3), dtype=np.cdouble) + #rho = Cube("../1.1-pp/Rho.cube") #just for struct + #rho_r_FFTbox[0,:,:,:] = rho.data + print('sum(rho_r)*(vol/nfft)',np.sum(rho_r_FFTbox)*(vol/nfft)) # debug - + rho_g_FFTbox = fft_r2g(rho_r_FFTbox) + rho_g, igvec = flat_FFTbox(rho_g_FFTbox) vh_r = get_vh(rho_g, igvec, bcell) + rho_r = rho_r_FFTbox + Eh = np.real(0.5*np.sum(vh_r*rho_r*(vol/nfft))) - vbh_r = vbh.data/2 # in Hartree + vbh_r = vbh.data/2 # from Ry to Hartree + Eh_ref = np.real(0.5*np.sum(vbh_r*rho_r*(vol/nfft))) diff = np.abs(vbh_r-vh_r) - print('Error = abs(vsub - vsub_aprox) in Hartree') + print(f'Eh_ref = {2*Eh_ref} in Ry') + print(f'Eh_ref-Eh = {2*(Eh_ref-Eh)} in Ry') + print(f'Eh_ref-Eh/Eh_ref = {(Eh_ref-Eh)/Eh_ref} ') + print('Error = abs(vh_qe - vh) in Hartree') print('max(Error)',np.max(diff)) print('min(Error)',np.min(diff)) print('avg(Error)',np.average(diff)) diff = np.abs(vbh_r/(vh_r)) - print('Ratio = abs(vbh_r/vh_r) in Hartree') + print('Ratio = abs(vbh_r/vh_r)') print('max(Error)',np.max(diff)) print('min(Error)',np.min(diff)) print('avg(Error)',np.average(diff)) diff --git a/WS2/LDA_SG15/4.0-hartree/density2Vh.py b/WS2/LDA_SG15/4.0-hartree/density2Vh.py new file mode 100644 index 0000000..fe9d877 --- /dev/null +++ b/WS2/LDA_SG15/4.0-hartree/density2Vh.py @@ -0,0 +1,172 @@ +import numpy as np +from fft_mod import put_FFTbox2, ifft_g2r, fft_r2g, flat_FFTbox +from pymatgen.io.cube import Cube +bohr2ang = 0.529177249 + + +def main(): + n1, n2, n3 = [18, 18, 135] + nfft = n1*n2*n3 + + # vh + vbh = Cube("../1.1-pp/Vh.cube") #just for struct + struct = vbh.structure + vol = vbh.volume/(bohr2ang**3) # in Bohr^3 + vol = 1449.7341 + cell = struct.lattice._matrix*(1/bohr2ang) # in Bohr + bcell = 2*np.pi*np.linalg.inv(cell).T # in Bohr^-1 + + # rho_r + rho_r_FFTbox = np.zeros((1,n1,n2,n3), dtype=np.cdouble) + rho = Cube("../1.1-pp/Rho.cube") #just for struct + rho_r_FFTbox[0,:,:,:] = rho.data + print('sum(rho_r)*(vol/nfft)',np.sum(rho_r_FFTbox)*(vol/nfft)) + # debug + rho_g_FFTbox = fft_r2g(rho_r_FFTbox) + rho_g, igvec = flat_FFTbox(rho_g_FFTbox) + vh_r = get_vh(rho_g, igvec, bcell) + rho_r = rho_r_FFTbox + Eh = np.real(0.5*np.sum(vh_r*rho_r*(vol/nfft))) + + vbh_r = vbh.data/2 # from Ry to Hartree + Eh_ref = np.real(0.5*np.sum(vbh_r*rho_r*(vol/nfft))) + + diff = np.abs(vbh_r-vh_r) + print(f'Eh_ref = {2*Eh_ref} in Ry') + print(f'Eh_ref-Eh = {2*(Eh_ref-Eh)} in Ry') + print(f'Eh_ref-Eh/Eh_ref = {(Eh_ref-Eh)/Eh_ref} ') + print('Error = abs(vsub - vsub_aprox) in Hartree') + print('max(Error)',np.max(diff)) + print('min(Error)',np.min(diff)) + print('avg(Error)',np.average(diff)) + diff = np.abs(vbh_r/(vh_r)) + print('Ratio = abs(vbh_r/vh_r)') + print('max(Error)',np.max(diff)) + print('min(Error)',np.min(diff)) + print('avg(Error)',np.average(diff)) + #print() + #print('diff = abs(vsub - vsub_aprox) in Hartree') + #print('max(diff)',np.max(diff)) + #print('min(diff)',np.min(diff)) + #print('avg(diff)',np.average(diff)) + #print(np.sort(ratio.flat)[-1000:]) + return + + +def get_vh(rho_g, igvec, bcell): + """ + --Input-- + rho_g : complex128(ng) + momentum space charge density in Hartree Atomic Unit + + gvec : int(ng,3) + momentum space charge density in Hartree Atomic Unit + + bcell_in: float64(3,3) + reciprocal lattice vectors (bohr^{-1}) + + --Output-- + vh : float(:,:,:) + + """ + vh_g = np.zeros_like(rho_g) + gvec = igvec.dot(bcell) + gvec_sq = np.einsum('ix,ix->i', gvec, gvec) # bohr^{-1} or Ry + for ig, g_sq in enumerate(gvec_sq): + if g_sq < 1e-12: + vh_g[ig] = 0 + else: + vh_g[ig] = 4*np.pi*rho_g[ig]/g_sq + + n1, n2, n3 = [18, 18, 135] + vh_g_FFTbox = put_FFTbox2(vh_g, igvec, [n1, n2, n3], noncolin=False) + vh_r = ifft_g2r(vh_g_FFTbox) + return np.real(vh_r[0,:,:,:]) + + + +def get_vxc(rho): + """ + --Input-- + rho : float(:,:,:) + real space charge density in Hartree Atomic Unit + + --Output-- + vxc : float(:,:,:) + + """ + inp = {} + n1, n2, n3 = rho.shape + rho_flatten = rho.flatten() + inp["rho"] = rho_flatten + + # Build functional + func_c = pylibxc.LibXCFunctional("LDA_C_PZ", "unpolarized") + func_x = pylibxc.LibXCFunctional("LDA_X", "unpolarized") + + # Compute exchange part + ret_x = func_x.compute(inp, do_fxc=True) + v2rho2_x = ret_x['v2rho2'][:,0] + vrho_x = ret_x['vrho'][:,0] + zk_x = ret_x['zk'][:,0] + + # Compute correlation part + ret_c = func_c.compute(inp, do_fxc=True) + v2rho2_c = ret_x['v2rho2'][:,0] + vrho_c = ret_c['vrho'][:,0] + zk_c = ret_c['zk'][:,0] + + vxc = vrho_x + vrho_c + fxc = v2rho2_x + v2rho2_c + + vxc_FFTbox = np.reshape(vxc,(n1,n2,n3)) + fxc_FFTbox = np.reshape(fxc,(n1,n2,n3)) + + return vxc_FFTbox, fxc_FFTbox + +def vxc_libxc(): + n1, n2, n3 = [18, 18, 135] + #gvec_rho = np.loadtxt('../2.1-wfn/rho_gvec.txt',dtype=np.int32) + #rho_pw2bgw = np.load('../2.1-wfn/rho_pw2bgw.npy') + #rho_FFTbox = put_FFTbox2(rho_pw2bgw, gvec_rho, [n1,n2,n3], noncolin=False) + #rho = ifft_g2r(rho_FFTbox) + #rho_pp = Cube("../1.1-pp/Rho.cube") + #rho = rho_pp.data + rho_pw2bgw = np.load("../2.1-wfn/rhor_pw2bgw.npy") + rho = rho_pw2bgw[0,:,:,:] + frac = 0.0001 + rho_f = rho * frac + rho_r = rho - rho_f + + # Create input + const = 1 + inp = {} + rho_flatten = rho.flatten() + rho_r_flatten = rho_r.flatten() + rho_f_flatten = rho_f.flatten() + inp["rho"] = const*rho_flatten + + # Build functional + func_c = pylibxc.LibXCFunctional("LDA_C_PZ", "unpolarized") + func_x = pylibxc.LibXCFunctional("LDA_X", "unpolarized") + + # Compute exchange part + ret_x = func_x.compute(inp, do_fxc=True) + v2rho2_x = ret_x['v2rho2'][:,0] + vrho_x = ret_x['vrho'][:,0] + zk_x = ret_x['zk'][:,0] + + # Compute correlation part + ret_c = func_c.compute(inp, do_fxc=True) + v2rho2_c = ret_x['v2rho2'][:,0] + vrho_c = ret_c['vrho'][:,0] + zk_c = ret_c['zk'][:,0] + + vxc = vrho_x + vrho_c + fxc = v2rho2_x + v2rho2_c + + vxc_FFTbox = np.reshape(vxc,(18,18,135)) + fxc_FFTbox = np.reshape(fxc,(18,18,135)) + return vxc_FFTbox, fxc_FFTbox + +main() diff --git a/WS2/LDA_SG15/4.0-hartree/fft_mod.py b/WS2/LDA_SG15/4.0-hartree/fft_mod.py index cb84341..284b64b 100644 --- a/WS2/LDA_SG15/4.0-hartree/fft_mod.py +++ b/WS2/LDA_SG15/4.0-hartree/fft_mod.py @@ -1,7 +1,54 @@ import os import time import numpy as np -from numpy.fft import fftn, ifftn, ifftshift +from numpy.fft import fftn, ifftn, fftshift, ifftshift + + +def flat_FFTbox(unkg_FFTbox): + """ + Generate a flattened 1d array 'unkg' and corresponding igvecs from unkg_FFTbox. + Assume unkg_FFTbox has its zero-frequency component in unkg_FFTbox[ispin,0,0,0]. + This functions is a kind of the inverse function of 'put_FFTbox'. + Note, because of the zero padding in unkg_FFTbox, unkg may have lots of zeros. + + --Input-- + unkg_FFTbox : cdouble(nspin,N1,N2,N3) + unkg in FFTbox + for colinear case, nspin = 1 + for non-colinear case, nspin = 2 + + --Output-- + unkg : cdouble(nspin*N1*N2*N3) + for non-colinear case, the first half is for spin-up, + the second half is for spin-dw. + gvecs_sym : int(N1*N2*N3,3) + array of miller index (h,k,l) + """ + nspin, n1, n2, n3 = unkg_FFTbox.shape + ng = n1*n2*n3 + h_max, k_max, l_max = n1//2, n2//2, n3//2 + + hs, ks, ls = np.unravel_index(np.array(range(ng)),(n1,n2,n3)) + + #generate zero starting positive index. it is just for reshaping + gvecs_pos = np.concatenate([[hs],[ks],[ls]], axis = 0).T + hkl2flat = ( n3 * (n2 * gvecs_pos[:, 0] + gvecs_pos[:, 1]) \ + + gvecs_pos[:, 2]) + hs[hs>=h_max] -= n1 + ks[ks>=k_max] -= n2 + ls[ls>=l_max] -= n3 + + #generate zero starting symmetric index. it will be returned + gvecs_sym = np.concatenate([[hs],[ks],[ls]], axis = 0).T + + unkg = np.zeros(nspin*ng,dtype=np.cdouble) + for ispin in range(nspin): + gstart = ispin*ng + gend = (ispin + 1)*ng + unkg[gstart:gend] = unkg_FFTbox[ispin,:,:,:].flat[hkl2flat] + + return unkg, gvecs_sym + def fft_r2g(unkr_FFTbox): """ fft from R space to G space @@ -22,7 +69,7 @@ def fft_r2g(unkr_FFTbox): for ispin in range(nspin): unkg_FFTbox[ispin,:,:,:] = fftn(unkr_FFTbox[ispin,:,:,:]) end = time.time() - print(f'\nIFFT time', end-start) + print(f'\nFFT time', end-start) return unkg_FFTbox @@ -33,6 +80,9 @@ def ifft_g2r(unkg_FFTbox): with an input array x[l], the output array y[k] satisfy y[k] = np.sum(x * np.exp(2j * np.pi * k * np.arange(n)/n)) / len(x) + Note unkg_FFTbox have the unk(G=0) component as the first element + in the unkg_FFTbox[ispin,0,0,0]. + ---Input--- unkg_FFTbox : cdouble(ispin,N1,N2,N3) @@ -76,6 +126,10 @@ def put_FFTbox2(unkg, gvecs, FFTgrid, noncolin, savedir='./', savefn=False): for colinear case, nspin = 1 for non-colinear case, nspin = 2 + unkg_FFTbox[ispin,0,0,0] = unk_ispin(G=0), + which is a zero frequency start array. + + """ #fn = savedir+f'./Unkg_FFTbox_k{ik}.npy' @@ -170,10 +224,3 @@ def put_FFTbox(ik, FFTgrid, savedir='./'): print(f'\nFFTBOX time', end-start) return unkg_FFTbox - -if __name__=='__main__': - ik = 1 - FFTgrid = [495, 495, 83] - unkg_FFTbox1 = put_FFTbox(ik, FFTgrid, savedir='../../data/') - unkg_FFTbox2 = put_FFTbox2(ik, FFTgrid, savedir='../../data/') - print(np.allclose(unkg_FFTbox2,unkg_FFTbox1)) diff --git a/WS2/LDA_SG15/4.0-hartree/fft_test.py b/WS2/LDA_SG15/4.0-hartree/fft_test.py new file mode 100644 index 0000000..c6dfd6e --- /dev/null +++ b/WS2/LDA_SG15/4.0-hartree/fft_test.py @@ -0,0 +1,165 @@ +import numpy as np +from fft_mod import put_FFTbox2, ifft_g2r, fft_r2g, flat_FFTbox +from pymatgen.io.cube import Cube +bohr2ang = 0.529177249 + + +def main(): + #rho_g = np.load("../2.1-wfn/rho_pw2bgw.npy") + #igvec = np.load("../2.1-wfn/rho_gvec.npy") + rho_g = np.load('../1-scf/rhog.npy') + igvec = np.load('../1-scf/igvec.npy') + n1, n2, n3 = [18, 18, 135] + nfft = n1*n2*n3 + + #rho = Cube("../1.1-pp/Rho.cube") #just for struct + #vbh = Cube("../1.1-pp/Vh.cube") #just for struct + #struct = vbh.structure + #vol = vbh.volume/(bohr2ang**3) # in Bohr^3 + #vol = 1449.7341 + #cell = struct.lattice._matrix*(1/bohr2ang) # in Bohr + #bcell = 2*np.pi*np.linalg.inv(cell).T # in Bohr^-1 + + # debug + rho_g_FFTbox = put_FFTbox2(rho_g, igvec, [n1, n2, n3], noncolin=False) # zero_init + rho_r = ifft_g2r(rho_g_FFTbox) + rho_g_test = fft_r2g(rho_r) # zero_init + rho_g_test_flatten, gvecs = flat_FFTbox(rho_g_test) + rho_g_FFTbox_2 = put_FFTbox2(rho_g_test_flatten, gvecs, [n1, n2, n3], noncolin=False) # zero_init + + + + + + # Check if f(g) --ifft--> g(r) --fft--> h(g) == f(g) + + diff = np.abs(rho_g_FFTbox_2 - rho_g_test) + print('Error = abs(rho_g_FFTbox - rho_g_test) ') + print('max(Error)',np.max(diff)) + print('min(Error)',np.min(diff)) + print('avg(Error)',np.average(diff)) + #diff = np.abs(vbh_r/(vh_r)) + #print('Ratio = abs(vbh_r/vh_r)') + #print('max(Error)',np.max(diff)) + #print('min(Error)',np.min(diff)) + #print('avg(Error)',np.average(diff)) + return + + +def get_vh(rho_g, igvec, bcell): + """ + --Input-- + rho_g : complex128(ng) + momentum space charge density in Hartree Atomic Unit + + gvec : int(ng,3) + momentum space charge density in Hartree Atomic Unit + + bcell_in: float64(3,3) + reciprocal lattice vectors (bohr^{-1}) + + --Output-- + vh : float(:,:,:) + + """ + vh_g = np.zeros_like(rho_g) + gvec = igvec.dot(bcell) + gvec_sq = np.einsum('ix,ix->i', gvec, gvec) # bohr^{-1} or Ry + for ig, g_sq in enumerate(gvec_sq): + if g_sq < 1e-12: + vh_g[ig] = 0 + else: + vh_g[ig] = 4*np.pi*rho_g[ig]/g_sq + + n1, n2, n3 = [18, 18, 135] + vh_g_FFTbox = put_FFTbox2(vh_g, igvec, [n1, n2, n3], noncolin=False) + vh_r = ifft_g2r(vh_g_FFTbox) + return np.real(vh_r[0,:,:,:]) + + + +def get_vxc(rho): + """ + --Input-- + rho : float(:,:,:) + real space charge density in Hartree Atomic Unit + + --Output-- + vxc : float(:,:,:) + + """ + inp = {} + n1, n2, n3 = rho.shape + rho_flatten = rho.flatten() + inp["rho"] = rho_flatten + + # Build functional + func_c = pylibxc.LibXCFunctional("LDA_C_PZ", "unpolarized") + func_x = pylibxc.LibXCFunctional("LDA_X", "unpolarized") + + # Compute exchange part + ret_x = func_x.compute(inp, do_fxc=True) + v2rho2_x = ret_x['v2rho2'][:,0] + vrho_x = ret_x['vrho'][:,0] + zk_x = ret_x['zk'][:,0] + + # Compute correlation part + ret_c = func_c.compute(inp, do_fxc=True) + v2rho2_c = ret_x['v2rho2'][:,0] + vrho_c = ret_c['vrho'][:,0] + zk_c = ret_c['zk'][:,0] + + vxc = vrho_x + vrho_c + fxc = v2rho2_x + v2rho2_c + + vxc_FFTbox = np.reshape(vxc,(n1,n2,n3)) + fxc_FFTbox = np.reshape(fxc,(n1,n2,n3)) + + return vxc_FFTbox, fxc_FFTbox + +def vxc_libxc(): + n1, n2, n3 = [18, 18, 135] + #gvec_rho = np.loadtxt('../2.1-wfn/rho_gvec.txt',dtype=np.int32) + #rho_pw2bgw = np.load('../2.1-wfn/rho_pw2bgw.npy') + #rho_FFTbox = put_FFTbox2(rho_pw2bgw, gvec_rho, [n1,n2,n3], noncolin=False) + #rho = ifft_g2r(rho_FFTbox) + #rho_pp = Cube("../1.1-pp/Rho.cube") + #rho = rho_pp.data + rho_pw2bgw = np.load("../2.1-wfn/rhor_pw2bgw.npy") + rho = rho_pw2bgw[0,:,:,:] + frac = 0.0001 + rho_f = rho * frac + rho_r = rho - rho_f + + # Create input + const = 1 + inp = {} + rho_flatten = rho.flatten() + rho_r_flatten = rho_r.flatten() + rho_f_flatten = rho_f.flatten() + inp["rho"] = const*rho_flatten + + # Build functional + func_c = pylibxc.LibXCFunctional("LDA_C_PZ", "unpolarized") + func_x = pylibxc.LibXCFunctional("LDA_X", "unpolarized") + + # Compute exchange part + ret_x = func_x.compute(inp, do_fxc=True) + v2rho2_x = ret_x['v2rho2'][:,0] + vrho_x = ret_x['vrho'][:,0] + zk_x = ret_x['zk'][:,0] + + # Compute correlation part + ret_c = func_c.compute(inp, do_fxc=True) + v2rho2_c = ret_x['v2rho2'][:,0] + vrho_c = ret_c['vrho'][:,0] + zk_c = ret_c['zk'][:,0] + + vxc = vrho_x + vrho_c + fxc = v2rho2_x + v2rho2_c + + vxc_FFTbox = np.reshape(vxc,(18,18,135)) + fxc_FFTbox = np.reshape(fxc,(18,18,135)) + return vxc_FFTbox, fxc_FFTbox + +main() diff --git a/WS2/LDA_SG15/4.0-hartree/fft_tutorial.py b/WS2/LDA_SG15/4.0-hartree/fft_tutorial.py new file mode 100644 index 0000000..0dc7df2 --- /dev/null +++ b/WS2/LDA_SG15/4.0-hartree/fft_tutorial.py @@ -0,0 +1,35 @@ +import numpy as np +from numpy.fft import fftn, ifftn, fftshift, ifftshift + + +sym_zero_cent = np.arange(-5, 6) +n1 = len(sym_zero_cent) +print('sym_zero_cent', sym_zero_cent) +sym_zero_init = ifftshift(sym_zero_cent) +print('sym_zero_init', sym_zero_init) + +pos_zero_init = np.arange(0,n1) +print('pos_zero_init', pos_zero_init) +pos_zero_cent = fftshift(pos_zero_init) +print('pos_zero_cent', pos_zero_cent) +shifted_pos_zero_init = (pos_zero_init - n1//2) # becomes sym zero_centered +print(np.allclose(shifted_pos_zero_init, sym_zero_cent)) +print(n1//2) +print('shifted_pos_zero_init', shifted_pos_zero_init) +print('ifftshift(shifted_pos_zero_init)', ifftshift(shifted_pos_zero_init)) +print('pos_zero_init', pos_zero_init) + +n1, n2, n3 = [4, 4, 4] +Nt = n1*n2*n3 +h_max, k_max, l_max = n1//2, n2//2, n3//2 +#generate zero initial positive index +hs, ls, ks = np.unravel_index(np.array(range(Nt)),(n1,n2,n3)) +#generate zero centered positive index +print(hs) +hs[hs>=h_max] -= n1 +print(hs) +hs = (hs-h_max) +ls = (ls-l_max) +ks = (ks-k_max) +hkls = np.concatenate([[hs],[ls],[ks]], axis = 0).T +