-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Woochang Kim
committed
Jul 20, 2023
1 parent
c833ee9
commit 73cb1e9
Showing
7 changed files
with
520 additions
and
31 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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() |
Oops, something went wrong.