green ftn 1

Makefile

@@ -9,9 +9,10 @@ default: main
 include arch.make
 
 # DMFT_OBJS = ed_operators.o
-DMFT_OBJS = fft.o frprmn.o ed_io.o ed_solver.o ed_config.o ed_hamiltonian.o ed_utils.o ed_basis.o ed_operators.o func.o
+DMFT_OBJS = fft.o frprmn.o ed_io.o ed_solver.o ed_config.o ed_hamiltonian.o ed_utils.o ed_basis.o ed_operators.o func.o \
+			ed_green.o
 
-MOD_OBJS = alloc.o memory.o memoryinfo.o sys.o parallel_params.o precision.o timer.o timestamp2.o ionew.o
+MOD_OBJS = alloc.o memory.o memoryinfo.o sys.o parallel_params.o precision.o timer.o timestamp2.o ionew.o 
 OBJS = bsd.o main.o $(DMFT_OBJS)
 
 COM_OBJS=$(OBJS) $(SYSOBJ)
@@ -28,7 +29,7 @@ $(FDF):
 # Dependencies
 main.o ed_config.o: $(FDF)
 
-main.o: parallel_params.o ed_config.o ed_hamiltonian.o ed_solver.o ionew.o
+main.o: parallel_params.o ed_config.o ed_hamiltonian.o ed_solver.o ionew.o ed_green.o
 alloc.o: parallel_params.o precision.o sys.o ionew.o
 ed_solver.o: ed_utils.o ed_hamiltonian.o ed_io.o
 ed_hamiltonian.o: ed_utils.o ed_basis.o ed_operators.o

ed_basis.f90

@@ -12,7 +12,6 @@ module ed_basis
     integer, pointer :: nbasis_up(:)
     integer, pointer :: nbasis_down(:)
     integer, pointer :: nbasis(:)       
-    integer :: nbasis_loc ! number of basis in the sector local to the node.
 
     ! For up,down basis, 4-bit integer is sufficient,
     ! because the maximum number of sites will not be more than 31.
@@ -29,6 +28,10 @@ module ed_basis
 
     subroutine ed_basis_init
         integer :: i, nd
+
+        if (node.eq.0) then
+            write(6,"(A)") "ed_basis: basis initiailization"
+        endif
 
         call re_alloc(nlocals,0,Nodes-1,name="nlocals",routine="ed_hamiltonian_init")
         call re_alloc(offsets,0,Nodes-1,name="offsets",routine="ed_hamiltonian_init")
@@ -57,8 +60,9 @@ subroutine ed_basis_init
         endif
     end subroutine ed_basis_init
 
-    subroutine prepare_basis_for_sector(isec)
+    subroutine prepare_basis_for_sector(isec,nbasis_loc)
         include 'mpif.h'
+        integer :: nbasis_loc ! number of basis in the sector local to the node.
         integer, intent(in) :: isec
         integer :: nam, i
         isector = isec
@@ -153,7 +157,8 @@ integer function get_basis_idx(basis)
     end function get_basis_idx
 
 #ifdef DEBUG
-    subroutine dump_basis
+    subroutine dump_basis(nbasis_loc)
+        integer :: nbasis_loc
         integer :: i,iunit,j
         integer(kind=kind_basis) :: basis
         character(len=100) :: fname

ed_green.f90

@@ -0,0 +1,28 @@
+module ed_green
+    use precision
+    use parallel_params
+    use ed_io
+    use ed_config
+
+    implicit none
+
+contains
+
+    subroutine ed_calc_green_ftn(nev_calc)
+        integer, intent(in) :: nev_calc
+
+        real(dp), allocatable :: eigval(:), pev(:)
+
+        integer :: i
+
+        call import_eigval(nev_calc,eigval,pev)
+
+        if (node.eq.0) then
+            do i=1,nev_calc
+                print *, i,eigval(i),pev(i)
+            enddo
+        endif
+
+    end subroutine ed_calc_green_ftn
+
+end module ed_green

ed_hamiltonian.f90

@@ -55,10 +55,8 @@ subroutine ed_set_band_structure
         enddo
 
         if (Node.eq.0) then
-            write(6,*) " ************************ "
-            write(6,*) "     IMPURITY LEVELS      "
-            write(6,*) " ************************ "
 
+            write(6,*) "ed_hamiltonian: initial impurity levels"
             do i = 1, Norb
                 write(6,'(i2,3x,e)') i, ef(i)
             enddo
@@ -189,8 +187,8 @@ subroutine multiply_H(nglobal,nloc,A,B)
     end subroutine multiply_H
 
 #ifdef DEBUG
-    subroutine dump_hamiltonian(isec)
-        integer :: ib,jb
+    subroutine dump_hamiltonian(isec,nbasis_loc)
+        integer :: ib,jb,nbasis_loc
         integer(kind=kind_basis) :: basisi, basisj
 
         real(dp), allocatable :: H(:,:), A(:), B(:)

ed_io.F90

@@ -0,0 +1,101 @@
+module ed_io
+    use precision
+    use parallel_params
+    use ed_config
+    use sys
+
+    implicit none
+
+    integer, parameter, private :: EIGVEC_IUNIT_BASE=100
+    integer, parameter, private :: EIGVAL_IUNIT_BASE=101
+contains
+
+    subroutine export_eigvec(isector,iev,node,nloc,nev,eigvec)
+        integer, intent(in) :: isector, iev, node, nloc, nev
+        real(dp), intent(in) :: eigvec(nloc)
+
+        integer :: i,iunit
+        character(len=150) :: fname
+
+        call eigvec_filename(isector,iev,node,fname)
+
+        iunit = EIGVEC_IUNIT_BASE + node
+#ifdef DATA_PLAINTEXT
+        open(unit=iunit,file=fname,status="replace")
+        write(iunit,"(F10.3)") (eigvec(i),i=1,nloc)
+        close(iunit)
+#else
+        open(unit=iunit,file=fname,status="replace",form="unformatted")
+        write(iunit) (eigvec(i),i=1,nloc)
+        close(iunit)
+#endif
+    end subroutine
+
+    subroutine export_eigval(nev_calc,eigval_all,pev,ind)
+        integer, intent(in) :: nev_calc, ind(nsector*nev)
+        real(dp), intent(in) :: eigval_all(nev*nsector)
+        real(dp), intent(in) :: pev(nev*nsector)
+
+        integer :: k
+
+#ifdef DATA_PLAINTEXT
+        open(unit=EIGVAL_IUNIT_BASE,file="eigenvalues.dat",form="formatted",&
+            status="replace")
+        write(EIGVAL_IUNIT_BASE,"(I4)") nev_calc
+        write(EIGVAL_IUNIT_BASE,"(I4,2F10.3)") (ind(k),eigval_all(ind(k)),pev(ind(k)),k=1,nev_calc)
+#else
+        open(unit=EIGVAL_IUNIT_BASE,file="eigenvalues.dat",form="unformatted",&
+            status="replace")
+        write(EIGVAL_IUNIT_BASE) nev_calc
+        write(EIGVAL_IUNIT_BASE) (ind(k),eigval_all(ind(k)),pev(ind(k)),k=1,nev_calc)
+#endif
+        close(EIGVAL_IUNIT_BASE)
+
+    end subroutine export_eigval
+
+    subroutine import_eigval(nev_calc,eigval,pev)
+        integer, intent(in) :: nev_calc
+        real(dp), allocatable, intent(out) :: eigval(:), pev(:)
+
+        integer :: i, iunit, idx, nev_calc_read
+        iunit = EIGVAL_IUNIT_BASE
+#ifdef DATA_PLAINTEXT
+        open(unit=iunit,file="eigenvalues.dat",status="old",form="formatted")
+#else
+        open(unit=iunit,file="eigenvalues.dat",status="old",form="unformatted")
+#endif
+
+#ifdef DATA_PLAINTEXT
+        read(iunit,*) nev_calc_read
+#else
+        read(iunit) nev_calc_read
+#endif
+        if (nev_calc_read.ne.nev_calc) then
+            if (node.eq.0) then
+                write(6,*) "nev_calc mismatch."
+            endif
+            call die
+            return
+        endif
+
+        allocate(eigval(nev_calc),pev(nev_calc))
+        do i=1,nev_calc
+#ifdef DATA_PLAINTEXT
+            read(iunit,*) idx, eigval(i), pev(i)
+#else
+            read(iunit) idx, eigval(i), pev(i)
+#endif
+        enddo
+        close(iunit)
+
+    end subroutine import_eigval
+
+    subroutine eigvec_filename(isector,iev,inode,fname)
+        integer, intent(in) :: isector, iev, inode
+        character(len=150), intent(out) :: fname
+
+        write(fname,"(A7,I2.2,A1,I3.3,A1,I4.4)") "eigvec-",isector,"-",iev,"-",inode
+
+    end subroutine eigvec_filename
+
+end module ed_io

ed_solver.f90

@@ -12,34 +12,31 @@ module ed_solver
 
     implicit none
 
-    real(dp), pointer :: eigval(:), eigval_all(:)
-    real(dp), pointer :: eigvec(:,:)
 
     public
 contains
 
-    subroutine ed_solver_init
-
-        call re_alloc(eigval,1,nev,name="eigval",routine="ed_solve")
-        call re_alloc(eigval_all,1,nev*nsector,name="eigval_all",routine="ed_solve")
-
-    end subroutine ed_solver_init
-
     subroutine ed_solve(iloop,nev_calc)
         integer, intent(in) :: iloop
         integer, intent(out) :: nev_calc
-        integer :: isector,i
-        integer :: ind(nsector*nev)
+        integer :: isector,i,nbasis_loc
+        real(dp), pointer :: eigvec(:,:)
+
+        real(dp) :: eigval(nev), eigval_all(nev*nsector)
         real(dp) :: pev(nev*nsector)
+        integer :: ind(nsector*nev)
 
         call timer('ed_solve',1)
 
         do isector=1,nsector
-            call prepare_basis_for_sector(isector)
-            call re_alloc(eigvec,1,nbasis_loc,1,nev,"eigvec","ed_solve", &
+            if (node.eq.0) then
+                write(*,*) "ed_solver: iloop = ",iloop," , isector = ",isector
+            endif
+            call prepare_basis_for_sector(isector,nbasis_loc)
+            call re_alloc(eigvec,1,nbasis_loc,1,nev,name="eigvec",routine="ed_solve", &
                          copy=.false.,shrink=.true.)
 
-            call diag(isector)
+            call diag(isector,nbasis_loc,eigval,eigvec)
 
             do i=1,nev
                 eigval_all((isector-1)*nev+i) = eigval(i)
@@ -70,16 +67,19 @@ subroutine ed_solve(iloop,nev_calc)
         endif
 
         call mpi_barrier(comm,ierr)
+        call de_alloc(eigvec,name="eigvec",routine="ed_solve")
         call timer('ed_solve',2)
         return
     end subroutine ed_solve
 
-    subroutine diag(isector)
+    subroutine diag(isector, nbasis_loc, eigval, eigvec)
 #ifdef DEBUG
         include 'debug.h'
 #endif
         include 'stat.h'
-        integer, intent(in) :: isector
+        integer, intent(in) :: isector, nbasis_loc
+        real(dp), intent(out) :: eigvec(nbasis_loc,nev)
+        real(dp), intent(out) :: eigval(nev)
 
         integer  maxnloc,maxnev,maxncv,ldv
         parameter  (maxnloc=3500000,maxnev=40,maxncv=60,ldv=maxnloc)

fft.f90

@@ -0,0 +1,48 @@
+        subroutine fft_1d(N,in_fft,out_fft)
+
+        implicit none
+
+        include "fftw3.f" 
+
+        integer*8:: N, plan
+        double precision:: in_fft(N)
+        double complex:: out_fft(N/2+1)
+
+        call dfftw_plan_dft_r2c_1d(plan,N,in_fft,out_fft,FFTW_ESTIMATE)
+        call dfftw_execute_dft_r2c(plan,in_fft,out_fft)   ! cautious! dft_r2c! not dft
+        call dfftw_destroy_plan(plan)
+
+        return
+        end 
+
+        subroutine cfft_1d(N,in_fft,out_fft)
+
+        implicit none
+
+        include "fftw3.f"
+
+        integer*8::N, plan
+        double complex:: in_fft(N), out_fft(N)
+
+        call dfftw_plan_dft_1d(plan,N,in_fft,out_fft,FFTW_FORWARD,FFTW_ESTIMATE)
+        call dfftw_execute_dft(plan,in_fft,out_fft)
+        call dfftw_destroy_plan(plan)
+
+        return
+        end
+
+        subroutine cfft_1d_bk(N,in_fft,out_fft)
+
+        implicit none
+
+        include "fftw3.f"
+
+        integer*8::N, plan
+        double complex:: in_fft(N), out_fft(N)
+
+        call dfftw_plan_dft_1d(plan,N,in_fft,out_fft,FFTW_BACKWARD,FFTW_ESTIMATE)
+        call dfftw_execute_dft(plan,in_fft,out_fft)
+        call dfftw_destroy_plan(plan)
+
+        return
+        end