      module bound2b
       use grids_module
       use constants
       implicit none
       real*8,allocatable,dimension(:,:) :: A
       private :: A

      contains

      subroutine bound()
      use grids_module
      use system
      use channels 
      implicit none
      real*8,dimension(1:np*alphanp%nchmax) :: wf
      real*8 :: lambda1, lambda2, lambda
      real*8 :: E1,E2,E
      real*8 :: sumdg,sum1,fz
      integer :: i

      if (.not. allocated(A)) allocate(A(1:np*alphanp%nchmax,1:np*alphanp%nchmax))


       fz = 1.d0 ! use for sign of wave function

       E1=-2.0/hbarc
!!!!secant method to find E1, iteration method to solve lambda
       call eigenvalue(E1,lambda1,wf)
       

       E2=0.9*E1
       call eigenvalue(E2,lambda2,wf)


       do while(abs(E1-E2)>1.E-8)
          E=E2-(E1-E2)*(lambda2-1.0d0)/(lambda1-lambda2)
          if (E>0 .or. abs(E)<1.e-3) then
          write(*,*) "no bound state "
          exit
          end if
          
          call eigenvalue(E,lambda,wf)
          write(*,*) "E=", E*hbarc, "lambda=",lambda
          E1=E2
          lambda1=lambda2
          E2=E
          lambda2=lambda
       end do

        write(*,*) "bound state energy is", E*hbarc
        
        
        


       sumdg=0.d0
       do i=1,np*alphanp%nchmax
       sumdg=sumdg+wf(i)*wf(i)
       end do

       write (*,*)
       write (*,*) ' norm of eigenvector from dgeev :',sumdg

!       normalize wavefunction
      sum1 = 0.d0
       do i=1,np
       sum1=sum1+wf(i)*wf(i)*pw(i)*pp(i)*pp(i) + wf(i+np)*wf(i+np)*pw(i)*pp(i)*pp(i)
       enddo

      sum1=1.d0/sqrt(sum1)

      do i=1,np*alphanp%nchmax
      wf(i)=wf(i)*sum1
      end do

c        sign of wave function

      if (wf(3).lt.0.d0) fz=-1.d0

      do i=1,np
      wf(i)=wf(i)*fz
      write(13,*) pp(i), wf(i)
      end do
      
      
c     calculate d-state probability 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccc计算D-wave的概率ccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      
      write(*,*) "d-state probability is"
      
      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccc计算<\phi|T+V|\phi>ccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


     

      end subroutine
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine eigenvalue(E,lambda,wf)
      use grids_module
      use system
      use channels
      implicit none
      real*8 :: E,lambda
      integer :: n,i
      integer :: info,lwork,loc(1)
      real*8, dimension(1:np*alphanp%nchmax) :: wr,wi,wf
      real*8, dimension(1:4*np*alphanp%nchmax) :: work
      real*8,dimension(1:np*alphanp%nchmax,1:np*alphanp%nchmax) :: vl,vr
   

      call Amat_bound(E)
! prepare calling the LAPACK routine DGEEV
      lwork=4*np*alphanp%nchmax !  the value is taken from the documentation of DGEEV
! the routine will destroy the original array A
! and return the real parts of the eigenvalues in WR and the imaginary parts in WI
! VR(:,i) contain the eigenstate (=wave function) of the right eigenvector of eigenvalue i
! VL(:,i) contain the eigenstate of the left eigenvector of eigenvalue i


      
      call DGEEV('N','V',NP*alphanp%nchmax,A,NP*alphanp%nchmax,WR,WI,VL,NP*alphanp%nchmax,VR,NP*alphanp%nchmax,WORK,LWORK,INFO)
                                                  

      if (info/=0) then
         write(*,*) "problem happened when calling DGEEV"
         stop
      end if

! locate maximal eigenvector
! here assume that this eigenvector is real
! then return this eigenvalue in ETA
! and the eigenvector in the global variable WF (which should be of dimension NP)

C     do n=1,nex-1
C      loc=MAXLOC(wr) ! return value is type of array, must defined with array type
C      wr(loc(1))=0
C     end do
       loc=MAXLOC(wr)


       lambda=wr(loc(1))
       wf=vr(:,loc(1))

      end subroutine


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine Amat_bound(E)
       use system
       use grids_module
       use precision
       use potnp
       implicit none
       real*8 :: g0,E
       integer :: ip,jp,nch 
    
       A=0.0d0
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccc计算A矩阵cccccccccccccccccccccccccccccccc      
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       


      end subroutine
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


      end module
