Skip to content
This repository has been archived by the owner on Sep 30, 2021. It is now read-only.

Commit

Permalink
Validated hamiltonian matrix element generation.
Browse files Browse the repository at this point in the history
  • Loading branch information
Young Woo Choi committed Apr 9, 2016
1 parent 0e13521 commit d2ac44f
Show file tree
Hide file tree
Showing 13 changed files with 2,129 additions and 107 deletions.
Binary file added .DS_Store
View file
Binary file not shown.
4 changes: 2 additions & 2 deletions arch.make
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,14 @@
#
SIESTA_ARCH = ifort-mpich
FC = mpif90
FFLAGS = -traceback -O3 -fast -no-ipo
FFLAGS = -traceback -fast -no-ipo -xSSE4.2
LIBS = -mkl=sequential ./lib/ARPACK/libarpack_OSX.a \
./lib/ARPACK/parpack_MPI-OSX.a

################################################################

RANLIB = ranlib
DEFS = -DMPI -DDEBUG
DEFS = -DMPI #-DDEBUG
#
.F.o:
$(FC) -c $(FFLAGS) $(DEFS) $<
Expand Down
18 changes: 18 additions & 0 deletions debug.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
!
!\SCCS Information: @(#)
! FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2
!
! %---------------------------------%
! | See debug.doc for documentation |
! %---------------------------------%

integer logfil, ndigit, mgetv0, &
msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,&
mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,&
mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd

common /debug/ &
logfil, ndigit, mgetv0, &
msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,&
mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,&
mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
14 changes: 12 additions & 2 deletions ed_basis.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,14 +154,24 @@ end function get_basis_idx

#ifdef DEBUG
subroutine dump_basis
integer :: i,iunit
integer :: i,iunit,j
integer(kind=kind_basis) :: basis
character(len=100) :: fname

call io_assign(iunit)
write(fname,"(A,I1)") "basis.node.",node
open(unit=iunit,file=fname,status="replace")
do i=1,nbasis_loc
write(iunit,"(I5,B)") i,ed_basis_get(offsets(node)+i)
write(iunit,"(I10,3x)",advance="no") i
basis = ed_basis_get(offsets(node)+i)
do j=1,2*Nsite
if (BTEST(basis,j-1)) then
write(iunit,"(I1)",advance="no") 1
else
write(iunit,"(I1)",advance="no") 0
endif
enddo
write(iunit,*)
enddo
close(iunit)
end subroutine dump_basis
Expand Down
4 changes: 2 additions & 2 deletions ed_config.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,15 +68,15 @@ subroutine ed_read_options
call re_alloc(vk,1,norb,1,nbath,routine='ed_read_options',name='vk')
if (fdf_block('DMFT.Baths', bfdf)) then
i = 1
do while((fdf_bline(bfdf, pline)) .and. (i .le. nbath))
do while( (i .le. nbath) .and. (fdf_bline(bfdf, pline)))
ek(Norb+i) = fdf_breals(pline,1)
i = i + 1
enddo
endif

if (fdf_block('DMFT.BathCouplings', bfdf)) then
i = 1
do while((fdf_bline(bfdf, pline)) .and. (i .le. nbath))
do while( i.le.nbath .and. (fdf_bline(bfdf, pline)))
do j=1,norb
vk(j,i) = fdf_breals(pline,j)
enddo
Expand Down
43 changes: 41 additions & 2 deletions ed_hamiltonian.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,18 +85,18 @@ subroutine find_hk(kx,ky,Hik)
end subroutine find_hk

subroutine multiply_H(nglobal,nloc,A,B)

real(dp), intent(in) :: A(nloc)
integer, intent(in) :: nglobal, nloc
real(dp), intent(out) :: B(nloc)

real(dp) :: A_all(nglobal)
real(dp), allocatable :: A_all(:)
integer :: i,j,k
integer :: ispin,jspin,iorb,jorb,ibath
integer(kind=kind_basis) :: basis_i, basis_j

real(dp) :: coeff_sum, coeff

allocate(A_all(nglobal))
call mpi_allgatherv(A,nloc,mpi_double_precision,A_all,&
nlocals,offsets,mpi_double_precision,comm,ierr)

Expand Down Expand Up @@ -184,5 +184,44 @@ subroutine multiply_H(nglobal,nloc,A,B)
enddo
enddo
enddo iloop

deallocate(A_all)
end subroutine multiply_H

#ifdef DEBUG
subroutine dump_hamiltonian(isec)
integer :: ib,jb
integer(kind=kind_basis) :: basisi, basisj

real(dp), allocatable :: H(:,:), A(:), B(:)
integer :: isec

isector = isec
allocate(H(nbasis(isector),nbasis(isector)))
allocate(A(nbasis(isector)),B(nbasis(isector)))

do ib=1,nbasis(isector)
basisi = ed_basis_get(ib)
A(:) = 0.0_dp
B(:) = 0.0_dp

A(ib) = 1.0_dp

call multiply_H(nbasis(isector),nbasis_loc,A,B)

do jb=1,nbasis(isector)
H(ib,jb) = B(jb)
enddo
enddo

open(unit=77,file="hamiltonian.dump",status="replace")
do ib=1,nbasis(isector)
do jb=1,nbasis(isector)
write(77,"(F8.3)", advance="no") H(ib,jb)
enddo
write(77,*)
enddo
close(77)
end subroutine dump_hamiltonian
#endif
end module ed_hamiltonian
16 changes: 0 additions & 16 deletions ed_operators.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,22 +220,6 @@ subroutine pair_exchange(basis_in,iorb,jorb,basis_out,coeff)
return
end subroutine pair_exchange

subroutine dump_basis(basis)
integer(kind=kind_basis) :: basis
integer :: isite,ispin,bitidx
do ispin=1,2
do isite=1,Nsite
bitidx = get_bitidx(isite,ispin)
if (BTEST(basis,bitidx)) then
write(*,"(I1)",advance="no") 1
else
write(*,"(I1)",advance="no") 0
endif
enddo
enddo
write(*,*)
end subroutine

subroutine number_op(basis_in,isite,ispin,basis_out,sgn)
integer(kind=kind_basis), intent(in) :: basis_in
integer, intent(in) :: isite,ispin
Expand Down
112 changes: 29 additions & 83 deletions ed_solver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,10 @@ module ed_solver
real(dp), pointer :: eigval(:), eigval_all(:)
real(dp), pointer :: eigvec(:,:)

integer :: ncv
public
contains

subroutine ed_solver_init
ncv = 2*nev

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")
Expand All @@ -29,7 +27,6 @@ end subroutine ed_solver_init
subroutine ed_solve(iloop)
integer, intent(in) :: iloop
integer :: isector,i

call timer('ed_solve',1)

do isector=1,nsector
Expand All @@ -39,90 +36,41 @@ subroutine ed_solve(iloop)

call diag(isector)

do i=1,nev
eigval_all((isector-1)*nev+i) = eigval(i)
enddo
enddo

call timer('ed_solve',2)
end subroutine ed_solve

! subroutine diag(isector)
! integer, intent(in) :: isector

! integer maxnloc,maxnev,maxncv,ldv
! parameter (maxnloc=3500000,maxnev=40,maxncv=60,ldv=maxnloc)
! real(dp) :: v(ldv,maxncv), workl(maxncv*(maxncv+8)), &
! workd(3*maxnloc), d(maxncv,2), resid(maxnloc), &
! ax(maxnloc)
! logical :: select(maxncv)
! character, parameter :: bmat = 'I'
! character(len=2), parameter :: which = 'SA'
! integer :: iparam(11), ipntr(11), lworkl, info, ido, nconv, i, j
! real(dp) :: sigma, tol
! real(dp) :: pdnorm2
! external :: pdnorm2, daxpy

! lworkl = ncv*(ncv+8)
! tol = 0.0

! iparam(1) = 1 ! ishfts
! iparam(3) = 500 ! maxitr
! iparam(7) = 1 ! mode

! info = 0
! ido = 0

! do
! call pdsaupd( comm, ido, bmat, nbasis_loc, which, nev, tol, resid, &
! ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info )
! if (ido .eq. -1 .or. ido .eq. 1) then
! call multiply_H(nbasis(isector),nbasis_loc,workd(ipntr(1)),workd(ipntr(2)))
! else
! exit
! endif
! enddo
! if ( info .lt. 0 ) then
! if ( node.eq. 0 ) then
! print *, ' Error with pdsaupd, info = ', info
! print *, iparam(5)
! call die
! endif
! else
! call pdseupd(comm,.true.,'All', select,d,v,ldv,sigma, &
! bmat, nbasis_loc, which, nev, tol, resid, ncv, v, ldv, &
! iparam, ipntr, workd, workl, lworkl, ierr )
! if ( ierr .ne. 0) then
! if ( node .eq. 0 ) then
! print *, ' Error with pdseupd, info = ', ierr
! call die
! endif
! else
! nconv = iparam(5)
! do j=1, nconv
! call multiply_H(nbasis(isector),nbasis_loc,v(1:ldv,j),ax)
! call daxpy(nbasis_loc, -d(j,1), v(1,j), 1, ax, 1)
! d(j,2) = pdnorm2( comm, nbasis_loc, ax, 1 )
! enddo
! call pdmout(comm, 6, nconv, 2, d, ncv, -6, &
! 'Ritz values and direct residuals')
! end if
! endif

! do i=1,nev
! eigvec(:,i) = v(:,i)
! eigval(i) = d(i,1)
! enddo
! end subroutine diag
subroutine diag(isector)
#ifdef DEBUG
include 'debug.h'
#endif
include 'stat.h'
integer, intent(in) :: isector

integer maxnloc,maxnev,maxncv,ldv
parameter (maxnloc=3500000,maxnev=40,maxncv=60,ldv=maxnloc)
real(dp) :: v(ldv,maxncv), workl(maxncv*(maxncv+8)), &
workd(3*maxnloc), d(maxncv,2), resid(maxnloc), &
ax(maxnloc)
logical :: select(maxncv)
character, parameter :: bmat = 'I'
character(len=2), parameter :: which = 'SA'
integer :: iparam(11), ipntr(11), lworkl, info, ido, nconv, i, j
real(dp) :: resid(nbasis_loc), d(ncv,2), workd(3*nbasis_loc), &
workl(ncv*(ncv+8)), sigma, ax(nbasis_loc), tol
logical :: select(ncv)
integer :: iparam(11), ipntr(11), lworkl, info, ido, nconv, i, j, ncv
real(dp) :: sigma, tol
real(dp) :: pdnorm2
external :: pdnorm2, daxpy

#ifdef DEBUG
ndigit = -3
logfil = 6
msaupd = 1
#endif DEBUG
ncv = nev*2

lworkl = ncv*(ncv+8)
tol = 0.0

Expand All @@ -134,26 +82,23 @@ subroutine diag(isector)
ido = 0

do
write(*,*) "pdsaupd", ido, info
call pdsaupd( comm, ido, bmat, nbasis_loc, which, nev, tol, resid, &
ncv, eigvec, nbasis_loc, iparam, ipntr, workd, workl, lworkl, info )
ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info )
if (ido .eq. -1 .or. ido .eq. 1) then
call multiply_H(nbasis(isector),nbasis_loc,workd(ipntr(1)),workd(ipntr(2)))
else
exit
endif
enddo

if ( info .lt. 0 ) then
if ( node.eq. 0 ) then
print *, ' Error with pdsaupd, info = ', info
print *, iparam
print *, iparam(5)
call die
endif
else
write(*,*) node, "pdseupd"
call pdseupd(comm,.true.,'All', select,d,eigvec,nbasis_loc,sigma, &
bmat, nbasis_loc, which, nev, tol, resid, ncv, eigvec, nbasis_loc, &
call pdseupd(comm,.true.,'All', select,d,v,ldv,sigma, &
bmat, nbasis_loc, which, nev, tol, resid, ncv, v, ldv, &
iparam, ipntr, workd, workl, lworkl, ierr )
if ( ierr .ne. 0) then
if ( node .eq. 0 ) then
Expand All @@ -163,8 +108,8 @@ subroutine diag(isector)
else
nconv = iparam(5)
do j=1, nconv
call multiply_H(nbasis(isector),nbasis_loc,eigvec(1:nbasis_loc,j),ax)
call daxpy(nbasis_loc, -d(j,1), eigvec(1,j), 1, ax, 1)
call multiply_H(nbasis(isector),nbasis_loc,v(1:ldv,j),ax)
call daxpy(nbasis_loc, -d(j,1), v(1,j), 1, ax, 1)
d(j,2) = pdnorm2( comm, nbasis_loc, ax, 1 )
enddo
call pdmout(comm, 6, nconv, 2, d, ncv, -6, &
Expand All @@ -173,6 +118,7 @@ subroutine diag(isector)
endif

do i=1,nev
eigvec(:,i) = v(:,i)
eigval(i) = d(i,1)
enddo
end subroutine diag
Expand Down
Binary file modified fdf/libfdf.a
View file
Binary file not shown.
Binary file modified libfdf.a
View file
Binary file not shown.
Loading

0 comments on commit d2ac44f

Please sign in to comment.