add siesta example

README.md

@@ -19,7 +19,6 @@ In each directory, we compute
 
 We only consider collinear calculation.
 
-
 For method 2.
 
     we first do pp.x and print Vxc(r) = Vtot(r) - Vbare(r) - Vhatree(r) in a FFTgrid and then

WS2/LDA_SG15/1-scf/scf.in

@@ -1,6 +1,6 @@
 &control
    prefix = 'WS2'
-   pseudo_dir = '/global/u2/w/wkim94/pseudo/sg15_oncv_upf_2020-02-06/'
+   pseudo_dir = '/global/cfs/cdirs/mp149/gsharing/wkim/codes/ONCVPseudoPack/SG15_v1.2/UPF/LDA'
    !pseudo_dir = '/global/u2/w/wkim94/pseudo/pseudodojo/nc-sr-04_pw_standard_upf/'
    !pseudo_dir = '/work2/08702/wkim94/frontera/pseudo/pseudodojo/nc-sr-04_pbe_standard_upf/'
    calculation = 'scf'
@@ -33,8 +33,8 @@ CELL_PARAMETERS {angstrom}
 1.575000000000  2.7279800219    0.00 
 0.000000000000  0.000000000000  25.0 
 ATOMIC_SPECIES 
-   W 183.84 'W_ONCV_PBE-1.2.upf'
-   S 78.96  'S_ONCV_PBE-1.2.upf'
+   W 183.84 'W_LDA_SR.SG15v1.2.UPF'
+   S 78.96  'S_LDA_SR.SG15v1.2.UPF'
 ATOMIC_POSITIONS crystal
  W   0.33333333333333 0.33333333333333 0.5000000000 
  S   0.66666666666667 0.66666666666667 0.4360000000 

WS2/LDA_SG15/2.1-wfn/compute_vxcdat_G.py

@@ -10,11 +10,15 @@ def main():
 
     gvec_rho = np.loadtxt('./rho_gvec.txt',dtype=np.int32)
     rho_pw2bgw = np.load('./rho_pw2bgw.npy')
+    print(np.sum(rho_pw2bgw))
     #n1, n2, n3 = vxc.data.shape
     n1, n2, n3 = [18, 18, 135]
     rho_g_FFTbox = put_FFTbox2(rho_pw2bgw, gvec_rho, [n1,n2,n3], noncolin=False)
     rho_r        = ifft_g2r(rho_g_FFTbox)[0,:,:,:]
-    print(rho_r/rho_pp_data)
+    #print(rho_r/rho_pp_data)
+    print(np.min(np.abs(rho_r/rho_pp_data)))
+    print(np.average(np.abs(rho_r/rho_pp_data)))
+    print(np.max(np.abs(rho_r/rho_pp_data)))
     #vxc_diag = np.zeros(len(ibnd_lst),dtype=np.complex128)
 
     #for ibnd in ibnd_lst:

WS2/LDA_SG15/2.1-wfn/convertG2R.py

@@ -0,0 +1,10 @@
+import numpy as np
+from fft_mod import put_FFTbox2, ifft_g2r
+
+data_g = np.load('./vxc_pw2bgw.npy')
+gvec   = np.loadtxt('./vxc_gvec.txt',dtype=np.int32)
+fn = 'vxc_pw2bgw_FFTbox'
+n1, n2, n3 = 18, 18, 135
+data_g_FFTbox = put_FFTbox2(data_g, gvec, [n1,n2,n3], noncolin=False)
+data_r_FFTbox = ifft_g2r(data_g_FFTbox)
+np.save(fn, data_r_FFTbox)

WS2/LDA_SG15/2.1-wfn/pw2bgw.in

@@ -12,15 +12,15 @@
    !wfng_dk3 = 0.0
    !rhototg_flag = .true.
    !rhototg_file = 'RHOtot'
-   rhog_flag = .true.
-   rhog_file = 'RHO'
+   !rhog_flag = .true.
+   !rhog_file = 'RHO'
    vxcg_flag = .true.
    vxcg_file = 'VXC'
-   vxc_flag = .true.
-   vxc_file = 'vxc.dat'
-   vxc_integral = 'r'
-   vxc_diag_nmin = 1
-   vxc_diag_nmax = 20
-   vxc_offdiag_nmin = 1
-   vxc_offdiag_nmax = 20
+   !vxc_flag = .true.
+   !vxc_file = 'vxc.dat'
+   !vxc_integral = 'r'
+   !vxc_diag_nmin = 1
+   !vxc_diag_nmax = 20
+   !vxc_offdiag_nmin = 1
+   !vxc_offdiag_nmax = 20
 /

WS2/LDA_SG15/2.1-wfn/zerocore/WS2.save/wfc_h5tonpy.py

@@ -0,0 +1,97 @@
+#!/usr/bin/env python
+import h5py
+from numpy import *
+import os
+from os.path import exists
+from scipy.spatial import KDTree
+bohr = 0.529177211 #Angstrom
+
+
+def main():
+    ibnd_lst = list(range(1, 20+1))
+    ik_lst = list(range(1,36+1))
+    savedir = './npy.save'
+    for ik in ik_lst:
+        print('Working on ik:',ik)
+        unkg_blst,  miller_ids = h5tonpy(ik, ibnd_lst)
+        save(savedir+f'/Gvecs_k{ik}.npy', miller_ids)
+        save(savedir+f'/Unkg_k{ik}.npy', unkg_blst)
+    return None
+
+
+def h5tonpy(ik, ibnd_lst):
+    wfc_fn = 'WS2.save/wfc'+str(ik)+'.hdf5'.format(ik)
+    unkg_set, miller_ids, npol = read_wfc(wfc_fn)
+
+    return unkg_set[ibnd_lst[0]-1:ibnd_lst[-1],:], miller_ids
+
+
+
+def read_wfc(fn):
+
+    #Read wfc file
+    print("Read ",fn)
+    f = h5py.File(fn)
+
+    gamma_only = f.attrs['gamma_only'] # bool
+    igwx = f.attrs['igwx'] # int
+    ik = f.attrs['ik'] # int k-point index
+    ispin = f.attrs['ispin'] # int
+    nbnd = f.attrs['nbnd'] # int
+    ngw = f.attrs['ngw'] # int
+    npol = f.attrs['npol'] # int
+    scale_factor = f.attrs['scale_factor'] #float
+    xk = f.attrs['xk'] # 3 components array in Bohr^-1?
+
+    miller_ids = array(f['MillerIndices']) # (3 x 3)
+#    savetxt('miller_ids.{0:03d}.txt'.format(ik),miller_ids)
+    print("\nMiller indices are saved in a text file")
+    rec_latt = zeros((3,3))
+    rec_latt[0,:] = array(f['MillerIndices'].attrs['bg1'])
+    rec_latt[1,:] = array(f['MillerIndices'].attrs['bg2'])
+    rec_latt[2,:] = array(f['MillerIndices'].attrs['bg3'])
+
+    #rec_latt = (2*pi/alat)*rec_latt
+
+
+    evc = array(f['evc'])
+    if npol == 2:
+        print('\nReading non-colinear spinor wavefunction')
+        unkg_set = read_evc_non_colinear(evc, igwx, nbnd)
+    else:
+        unkg_set = read_evc(evc, igwx, nbnd)
+#    savetxt('unkg.{0:03d}.txt'.format(ik),unkg_set)
+#    print("\nWavefunctions are saved in a text file")
+    return unkg_set, miller_ids, npol
+
+
+def read_evc(evc, igwx, nbnd):
+
+    print('converting the wavefunction coefficents to numpy format')
+    psi_k_set = zeros((nbnd,igwx), dtype=complex64)
+    for n, row in enumerate(evc):
+        psi = add(row[::2], 1j*row[1::2])
+        psi_k_set[n,:] = psi
+
+    print('converting a wavefunction file is done')
+
+    return psi_k_set
+
+
+def read_evc_non_colinear(evc, igwx, nbnd):
+    #if npol == 2 the len of array will be doulbed
+
+    print('converting the wavefunction coefficents to numpy format')
+    psi_k_set = zeros((nbnd,igwx*2), dtype=complex64)
+    for n, row in enumerate(evc):
+        psi = add(row[::2], 1j*row[1::2])
+        psi_k_set[n,:] = psi
+
+    print('converting a wavefunction file is done')
+
+    return psi_k_set
+
+
+
+if __name__=='__main__':
+    main()