Added an option to output real part of WFs with phases.

This is available only for Gamma-point calculations.

tools/PostProcessing/local_module.f90

@@ -37,6 +37,7 @@ module local
   integer :: flag_proc_band_str
   logical :: flag_only_charge, flag_by_kpoint, flag_wf_range, flag_proc_range, flag_procwf_range_Ef
   logical :: flag_total_iDOS, flag_write_forces, flag_write_spin_moments, flag_l_resolved, flag_lm_resolved
+  logical :: flag_outputWF_real
   character(len=80) :: charge_stub
 
   integer :: i_job ! Job type

tools/PostProcessing/process_module.f90

@@ -96,7 +96,7 @@ contains
     use local, ONLY: nkp, efermi, current, nptsx, nptsy, nptsz, eigenvalues, flag_by_kpoint, &
          n_bands_total, band_active_kp, flag_proc_range, band_full_to_active, evec_coeff,&
          E_procwf_min, E_procwf_max, flag_procwf_range_Ef, band_proc_no, n_bands_process, &
-         grid_x, grid_y, grid_z, wtk
+         grid_x, grid_y, grid_z, wtk, flag_outputWF_real
     use output, ONLY: write_cube
     use global_module, only : nspin
     use read, ONLY: read_eigenvalues, read_psi_coeffs
@@ -104,6 +104,7 @@ contains
     use global_module, ONLY: ni_in_cell, atom_coord, species_glob
     use species_module, ONLY: nsf_species, n_species
     use angular_coeff_routines, ONLY: set_prefac_real, set_fact, set_prefac
+    use GenComms, ONLY: cq_abort
 
     implicit none
 
@@ -129,6 +130,10 @@ contains
     allocate(current(nptsx,nptsy,nptsz))
     allocate(psi(nptsx,nptsy,nptsz))
     max_band_integral_deviation = zero
+
+    if (flag_outputWF_real .and. (nkp.ne.1)) &
+       call cq_abort("OutputWF_real is available only for Gamma-point calculations.")
+
     if(flag_proc_range) then
        Emin = E_procwf_min
        Emax = E_procwf_max
@@ -154,7 +159,11 @@ contains
                       psi = zero
                       current = zero
                       call pao_to_grid(band_full_to_active(band), kp, ispin, psi)
-                      current = psi*conjg(psi)
+                      if (.not.flag_outputWF_real) then
+                         current = psi*conjg(psi)   ! band density
+                      else
+                         current = real(psi)        ! only real-part of WF
+                      endif
                       write(ci,'("Band",I0.6,"den_kp",I0.3,"S",I0.1)') band, kp, ispin
                       call write_cube(current,ci)
                       integral = gpv*sum(current)
@@ -183,7 +192,11 @@ contains
                       if(integral_deviation>band_integral_tol) &
                            write(*,fmt='(4x,"Integral of band ",i5," at kp ",i5," is ",f17.10)') &
                            band,kp,integral
-                      current = current + psi*conjg(psi)*wtk(kp)
+                      if (.not.flag_outputWF_real) then
+                         current = current + psi*conjg(psi)*wtk(kp)   ! band density
+                      else
+                         current = current + real(psi)*wtk(kp)        ! only real-part of WF
+                      endif
                       idum1 = 1
                    end if
                 end do ! kp
@@ -218,7 +231,11 @@ contains
                       psi = zero
                       current = zero
                       call pao_to_grid(band_full_to_active(band_proc_no(band)), kp, ispin, psi)
-                      current = psi*conjg(psi)
+                      if (.not.flag_outputWF_real) then
+                         current = psi*conjg(psi)   ! band density
+                      else
+                         current = real(psi)        ! only real-part of WF
+                      endif
                       write(ci,'("Band",I0.6,"den_kp",I0.3,"S",I0.1)') band_proc_no(band), kp, ispin
                       call write_cube(current,ci)
                       integral = gpv*sum(current)
@@ -245,7 +262,11 @@ contains
                       if(integral_deviation>band_integral_tol) &
                            write(*,fmt='(4x,"Integral of band ",i5," at kp ",i5," is ",f17.10)') &
                            band,kp,integral
-                      current = current + psi*conjg(psi)*wtk(kp)
+                      if (.not.flag_outputWF_real) then
+                         current = current + psi*conjg(psi)*wtk(kp)   ! band density
+                      else
+                         current = current + real(psi)*wtk(kp)        ! only real-part of WF
+                      endif
                    end if
                 end do ! kp
                 write(ci,'("Band",I0.6,"den_totS",I0.1)') band_proc_no(band), ispin

tools/PostProcessing/read_module.f90

@@ -255,6 +255,9 @@ contains
        stop
     end if
     flag_by_kpoint = fdf_boolean('Process.outputWF_by_kpoint',.false.)
+    ! if output only the real part of WFs (for Gamma-point only)
+    flag_outputWF_real = .false.
+    if (leqi(job,'ban')) flag_outputWF_real = fdf_boolean('Process.outputWF_real',.false.)
     ! DOS
     ! Add flag for window relative to Fermi level
     E_DOS_min = fdf_double('Process.min_DOS_E',E_wf_min)