Feature: Ground state projection tool for RT-TDDFT

tools/plot-tools/README.md

@@ -147,27 +147,3 @@ Then, the following command will output parsed partial DOS to directory `PDOS_FI
 abacus-plot -d -o
 ```
 
-### Dipole and Absorption
-
-```python
-
-from abacus_plot.dipole import Dipole
-from abacus_plot.dipole import Absorption
-import matplotlib.pyplot as plt
-
-dipolefile = './SPIN1_DIPOLE'
-dipole = Dipole(dipolefile, dt=0.0024)
-Efile=[["efield_0.dat"],["efield_1.dat"],["efield_2.dat"]]
-#Efile is a 2D list, Efile[i][j] is the jth Efield data file in ith direction
-Abs = Absorption(dipolefile, Efile, dt=0.0024)
-
-fig1, ax1 = plt.subplots()
-fig1, ax1 = dipole.plot_dipole(fig1, ax1)
-fig1.savefig('dipole.png')
-
-fig2, ax2 = plt.subplots()
-x_range = [0, 20]
-fig2, ax2 = Abs.plot_abs(
-    fig2, ax2, x_range=x_range, unit='eV')
-fig2.savefig('abs.png')
-```

tools/rt-tddft-tools/examples/ground-state-projection-Si/INPUT

@@ -0,0 +1,37 @@
+INPUT_PARAMETERS
+basis_type      lcao
+
+ecutwfc         100
+scf_nmax        100
+scf_thr         1e-6
+
+calculation     md
+esolver_type    tddft
+md_type         nve
+md_nstep        10000
+md_dt           0.002
+md_tfirst       0
+
+td_vext         1
+td_vext_dire    3
+td_stype        1
+td_ttype        0
+td_tstart               1
+td_tend                 5000
+td_lcut1                0.01
+td_lcut2                0.99
+td_gauss_freq           0.96
+td_gauss_phase          0.0
+td_gauss_sigma          0.5
+td_gauss_t0              1200
+td_gauss_amp            0.1
+
+out_dipole      1
+out_efield      1
+out_current     1
+out_current_k   1
+
+out_wfc_lcao    1
+out_mat_hs      1
+out_app_flag    0
+out_interval    25

tools/rt-tddft-tools/examples/ground-state-projection-Si/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+8 8 8 0 0 0

tools/rt-tddft-tools/examples/ground-state-projection-Si/STRU

@@ -0,0 +1,22 @@
+ATOMIC_SPECIES
+Si  28.085  Si_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+Si_gga_8au_100Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+10.2  // add lattice constant
+
+LATTICE_VECTORS
+0.5 0.5 0.0
+0.5 0.0 0.5
+0.0 0.5 0.5
+
+ATOMIC_POSITIONS
+Cartesian //Cartesian or Direct coordinate.
+
+Si      // Element type
+0.0     // magnetism
+2       // number of atoms
+0.00 0.00 0.00 m 0 0 0
+0.25 0.25 0.25 m 0 0 0

tools/rt-tddft-tools/examples/ground-state-projection-Si/projection.py

@@ -0,0 +1,111 @@
+import numpy as np
+import matplotlib.pyplot as plt
+class Projection:
+    def __init__(self, stepref, klist, steps, fdir='./OUT.ABACUS', wfc_dir='', s_dir=''):
+        self.stepref = stepref
+        self.klist = klist
+        self.steps = steps
+        self.wfc_dir = fdir + wfc_dir
+        self.s_dir = fdir + s_dir
+        wfc_ref, Ocp_ref = self.read_wfc(klist[0]+1, stepref+1, dir=self.wfc_dir)
+        self.nband = len(Ocp_ref)
+        self.nlocal = len(wfc_ref[0])
+
+    def read_wfc(self, kpoint, nstep, dir='.'):
+        file = dir+'/WFC_NAO_K'+str(kpoint)+'_ION'+str(nstep)+'.txt'
+        read_flag=0
+        with open(file,'r') as f:
+            for line in f:
+                if('bands' in line):
+                    bands=int(line.split()[0])
+                elif('orbitals' in line):
+                    orbitals=int(line.split()[0])
+                    wavef0=np.zeros([bands,orbitals],dtype=complex)
+                    Ocp=np.zeros(bands,dtype=float)
+                elif('(band)' in line):
+                    read_flag=0
+                    i=int(line.split()[0])-1
+                    j=0
+                elif('Occupations' in line):
+                    Ocp[i]=float(line.split()[0])
+                    read_flag=1
+                    continue
+                elif(read_flag==1):
+                    tmp=line.split()
+                    for l in range(len(tmp)//2):
+                        wavef0[i][j]=complex(float(tmp[2*l]),float(tmp[2*l+1]))
+                        j+=1
+        return wavef0,Ocp
+
+    def CSC(self, Cd,Sk,C):
+        Cdag=np.conjugate(Cd).transpose()
+        return np.matmul(Cdag,np.matmul(Sk,C))
+    def CSC2(self, Cd,Sk,C):
+        CSC1=self.CSC(Cd,Sk,C)
+        return CSC1*CSC1.conjugate()
+
+    def S_read(self, kpoint,nstep,dir='.'):
+        file = dir+'/'+str(nstep)+'_data-'+str(kpoint)+'-S'
+        count=0
+        fir=1
+        with open(file, 'r') as f1:
+            for line in f1:
+                tmp=line.split()
+                if(fir==1):
+                    dim=int(tmp[0])
+                    i=0
+                    Sk=np.zeros([dim,dim],dtype=complex)
+                    fir=0
+                    for j in range(dim):
+                        c_s=eval(tmp[j+1])
+                        Sk[i][j]=complex(c_s[0],c_s[1])
+                else:
+                    for j in range(dim-i):
+                        c_s=eval(tmp[j])
+                        Sk[i][i+j]=complex(c_s[0],c_s[1])
+                i+=1
+                if(i==0):
+                    break 
+        for i in range(dim):
+            for j in range(i):
+                Sk[i][j]=Sk[j][i].conjugate()
+        return Sk
+
+    def cal_Pnm(self, ik, nstep, wfc_ref):
+        wfc_t,Ocp_t=self.read_wfc(ik+1, nstep+1, dir=self.wfc_dir)
+        S_t=self.S_read(ik,nstep,dir=self.s_dir)
+        Pnm=np.zeros([self.nband,self.nband],dtype=float)
+        for a in range(self.nband):
+            for b in range(self.nband):
+                Pnm[a][b] = self.CSC2(wfc_ref[a],S_t,wfc_t[b]).real
+        return Pnm, Ocp_t
+
+    def cal_On_single(self, ik, nstep, wfc_ref):
+        Pnm, Ocp_t = self.cal_Pnm(ik, nstep, wfc_ref)
+        On = np.zeros(self.nband,dtype=float)
+        On = Ocp_t @ Pnm
+        return On
+
+    def save_On_all(self):
+        for ik in self.klist:
+            wfc_ref,Ocp_ref=self.read_wfc(ik+1, self.stepref+1, dir=self.wfc_dir)
+            On_tot = np.zeros([len(steps),self.nband],dtype=float)
+            for n, nstep in enumerate(steps):
+                On_tot[n] = self.cal_On_single(ik, nstep, wfc_ref)
+            np.savetxt('On_'+str(ik)+'.dat', On_tot)
+
+
+if __name__ == "__main__":
+    #the kpoints you need, check kpoints file to get the index, for this example, 0 means gamma point
+    klist=[0, 1]
+    #the steps you need, check STRU_MD file to get the index
+    start_step = 0
+    end_step = 10005
+    out_interval = 25
+    steps = range(start_step, end_step, out_interval)
+    #the ground state step
+    stepref=0
+    #Ground state projection
+    pro=Projection(0, klist, steps)
+    #save the Occupation infos to file
+    pro.save_On_all()