gpoda1ds.f90

Ground state of the time-independent Gross-Pitaevskii equation

!**********************************************************************
!
! Calculates the ground state of the time-independent 
! Gross-Pitaevskii equation using the optimal damping algorithm
!
! Spectral method in 1D
!
! Claude M. Dion and Eric Cances
!
! Version 1.0 (January 2007)
! 
!**********************************************************************

module criteria
  real(kind(1.d0)) :: critODA,critCG
end module criteria

module system1D
  integer :: sf
  real(kind(1.d0)) :: lambda
  real(kind(1.d0)), dimension(:,:), allocatable :: H0
  logical :: symmetric
end module system1D

module basis_set1D
  real(kind(1.d0)), dimension(:), allocatable :: expo
  real(kind(1.d0)), dimension(:,:), allocatable :: phi
end module basis_set1D


program GPODA1Ds
  !
  ! main program
  !
  use basis_set1D
  use system1D

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer :: n,Ng

  real(dp), dimension(:), allocatable :: C,psi2in,psi2out

  integer :: itMax,iter
  real(dp) :: mu,error,Eopt,slope,H0Cin,HinCin
  
  logical :: guess_from_file,output_grid,converged

  write(*,*) "GPODA1Ds"
  write(*,*)

  !
  ! Initialization
  !
  
  call read_params(n,Ng,itMax,guess_from_file,output_grid)
 
  write(*,*) "Initialization"

  allocate(phi(0:n,Ng),expo(Ng),H0(0:n,0:n))
  call initH(Ng,n)

  allocate(psi2in(Ng),psi2out(Ng),C(0:n))

  write(*,*)
  if(guess_from_file) then 
     ! read initial guess from file guess1Ds.data
     write(*,*) "  Read initial guess from file guess1Ds.data"
     call read_coeff(n,"guess1Ds.data",C)
     call basis2grid(n,Ng,C,psi2in)
     psi2in = psi2in**2
  else
     ! use ground state of H_0 as guess 
     write(*,*) "  Compute the ground state of H_0"
     C = 0.d0
     C(0) = 1.d0
     psi2in = 0.d0
     call ground_state(n,Ng,psi2in,mu,C,psi2in)
  endif
  write(*,*)

  ! compute initial values of H_0 C and H(psi) C
  call evalH0(n,Ng,C,C,H0Cin)
  call evalH(n,Ng,psi2in,C,C,HinCin)
  
  !
  ! Iterations of the ODA
  !
  do iter=1,ItMax
     write(*,*)
     write(*,*) "Iteration ",iter
     write(*,*) 
     write(*,*) "  Compute the ground state of H(psi_in)"
     call ground_state(n,Ng,psi2in,mu,C,psi2out)
     write(*,*) "  Optimal damping"
     call oda(n,Ng,C,psi2in,psi2out,H0Cin,HinCin,psi2in,Eopt,slope, &
          & H0Cin,HinCin,converged)
     write(*,*) 
     if(converged) exit
  end do

  if(converged) then
     write(*,*) 
     write(*,*) 'Convergence achieved in',iter,' iterations'
     write(*,*) ' --> mu =',mu
     write(*,*) ' --> E = ',0.5d0*(H0Cin+HinCin)
     write(*,*)
     
     write(*,*)
     write(*,*) "Checking self-consistency"
     write(*,*)         
     call ground_state(n,Ng,psi2in,mu,C,psi2out)
     call convergence_test_norm(Ng,psi2in,psi2out,converged,error)
     write(*,89) error,Ng,error/real(Ng,dp)
89   format("   l1 norm of psi2out-psi2in: ",e14.6," for ",i7," grid points"/ &
          & "   (",e12.6," per grid point)"/)
  else
     write(*,*) 'Not converged after ',itMax,' iterations'
  endif
     
  ! convention: lowest coefficient is positive
  if (C(0) < 0.d0) C = -C
  
  ! write out coefficients of the ground state
  call write_coeff(n,C,"gs1Ds.data")
  
  ! write out order parameter on a grid
  if(output_grid) call write_grid(n,C,"gs1Ds_grid.data")
  
  deallocate(H0,phi,expo,psi2in,psi2out,C)
  
end program GPODA1Ds


!**********************************************************************
!
! Input/output routines
!
!**********************************************************************

subroutine read_params(n,Ng,itMax,guess_from_file,output_grid)
  !  
  ! Reads parameters from the file params1Ds.in
  !
  use criteria
  use system1D

  implicit none
  integer, parameter :: dp=kind(1.d0)

  integer, intent(out) :: n,Ng
  integer, intent(out) :: itMax
  logical, intent(out) :: guess_from_file,output_grid

  namelist /params1Ds/ &
       & lambda, &          ! non-linear factor
       & n, &               ! highest basis function
       & symmetric, &       ! symmetric?
       & critODA, &         ! convergence criterion for the ODA
       & itMax, &           ! maximum number of iterations of the ODA
       & guess_from_file, & ! read initial guess from file guess1Ds.data?
       & output_grid        ! ouput order parameter on a grid?

  ! read parameters from file
  open(unit=8,file="params1Ds.in",status='old')
  read(8,params1Ds)
  close(8)
  
  ! check validity of some parameters
  if(n < 1) then
     write(*,*) "Invalid input parameter n =",n
     stop
  elseif(n > 91) then
     write(*,*) "Too many basis functions n =",n
     write(*,*) "Maximum value is 91"
     stop
  end if

  if(critODA <= 0.d0) then
     write(*,*) "Invalid criterion critODA =",critODA
     stop
  end if

  ! calculate constant to reconstruct HO number from index
  if(symmetric) then
     sf = 2
     n = n/2
  else
     sf = 1
  end if

  ! calculate number of spatial grid points
  Ng = 2*sf*n+1

  ! print out parameters 
  write(*,90) lambda,n+1,Ng

90 format("Parameters:"/,/," Nonlinearity =      ",e16.8,/,/,&
        & " Number of basis functions: ",i4/&
        & " Number of grid points:     ",i4)

  if(symmetric) then
     write(*,*) "Symmetric"
  end if
        
  write(*,*)

  return
end subroutine read_params


subroutine read_coeff(n,infile,c)
  !
  ! Read initial basis set coefficients
  !
  use system1D
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n
  character(len=*), intent(in) :: infile
  real(dp), dimension(0:n), intent(out) :: c

  integer :: i
  real(dp) :: norm,tmp

  ! Initialize coefficients
  c = 0.d0
  
  ! Read coefficients
  open(13,file=infile,status='old')
  reading:  do
     read(13,*,end=66) i,tmp
     
     if(symmetric) then
        if(2*(i/2) /= i) then
           write(*,*) 'Component ',i,' of input file ignored'
           cycle reading
        end if
        i = i/sf
     end if
     if(i > n .or. i < 0) then
        write(*,*) 'Component ',i,' of input file ignored'
        cycle reading
     end if
     
     if(c(i) /= 0.d0) then
        write(*,*) 'Warning!  Component ',i,' defined more than once,', &
             & ' last value kept'
     end if
     
     c(i) = tmp
  end do reading
  
66 close(13)

  ! Renormalize coefficients's
  norm = sum(c**2)
  c = c/sqrt(norm)

  return
end subroutine read_coeff


subroutine write_coeff(n,c,outfile)
  !
  ! Write basis set coefficients
  !
  use system1D
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n
  real(dp), dimension(0:n), intent(in) :: c
  character(len=*), intent(in) :: outfile

  integer :: i,ii

  open(unit=12,file=trim(outfile))

  do i=0,n
     ii = sf*i
     write(12,*) ii,C(i)
  enddo

  close(12)

  return
end subroutine write_coeff


subroutine write_grid(n,C,outfile)
  !
  ! Write order parameter on a grid
  !
  use system1D
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n
  real(dp), dimension(0:n), intent(in) :: C
  character(len=*), intent(in) :: outfile

  integer :: ng,i,ii
  real(dp) :: xmin,xmax,dx
  real(dp), dimension(:), allocatable :: x,phit,psi
  real(dp), dimension(:,:), allocatable :: phi

  namelist /grid1D/ &
       & ng, &   ! number of grid points
       & xmin, & ! first point 
       & xmax    ! last point 

  ! read grid parameters
  open(unit=8,file="params1Ds.in",status='old')
  read(8,grid1D)
  close(8)

  ! construct grid
  allocate(x(ng))

  dx = (xmax-xmin)/real(ng-1,dp)
  do ii = 1,ng
     x(ii) = xmin + real(ii-1,dp)*dx
  end do

  ! transform x coordinate
  allocate(phi(0:n,ng),phit(ng))
  phi = 0.d0

  do i = 0,n
     ii = sf*i
     call harmonic(ii,ng,x,phit)
     phi(i,:) = phit
  end do
  deallocate(phit)
  
  allocate(psi(ng))
  psi = 0.d0
  do ii = 1,ng
     do i = 0,n
        psi(ii) = psi(ii)+phi(i,ii)*C(i)
     end do
  end do
  deallocate(phi)

  ! write order parameter
  open(unit=12,file=trim(outfile))
  
  do ii = 1,ng
     write(12,88) x(ii),psi(ii)
88   format(e14.6,1x,e16.8)
  end do

  close(12)

  deallocate(x,psi)

  return
end subroutine write_grid


!**********************************************************************
!
! Minimization routines
!
!**********************************************************************

subroutine ground_state(n,Ng,psi2,eigenval,C,psi2out)
  !
  ! Calculate the ground state of H(psi) 
  ! by direct diagonalization
  !
  use criteria

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n,Ng
  real(dp), dimension(Ng), intent(in) :: psi2
  real(dp), intent(out) :: eigenval
  real(dp), dimension(0:n), intent(out) :: C
  real(dp), dimension(Ng), intent(out) :: psi2out
  
  real(dp), dimension(0:n,0:n) :: H
  real(dp), dimension(Ng) :: psi

  ! construct Hamiltonian matrix
  call makeH(n,Ng,psi2,H)  

  ! find eigenfunction with smallest eigenvalue
  call eigen(n+1,H,eigenval,C)

  write(*,*) "   --> mu =",eigenval

  call basis2grid(n,Ng,C,psi)
  psi2out = psi**2

  return
end subroutine ground_state


subroutine oda(n,Ng,Cout,psi2in,psi2out,H0Cin,HinCin,psi2,Eopt,slope, &
     & H0C,HC,converged)
  !
  ! Optimal Damping Algorithm
  !
  use criteria
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n,Ng
  real(dp), dimension(0:n), intent(in) :: Cout
  real(dp), dimension(Ng), intent(in) :: psi2in,psi2out
  real(dp), intent(in) :: H0Cin,HinCin
  real(dp), dimension(Ng), intent(out) :: psi2
  real(dp), intent(out) :: H0C,HC
  logical, intent(out) :: converged

  real(dp) :: H0Cout,HoutCout,HinCout
  real(dp) :: energyIn,slope,curv,alpha,Eopt
  
  call evalH0(n,Ng,Cout,Cout,H0Cout)
  call evalH(n,Ng,psi2out,Cout,Cout,HoutCout)
  call evalH(n,Ng,psi2in,Cout,Cout,HinCout)
  
  energyIn = 0.5d0*(H0Cin+HinCin)
  slope = HinCout - HinCin
  curv = HinCin + HoutCout - 2.d0*HinCout+H0Cout-H0Cin
  
  alpha = -slope/curv
  if(alpha >= 1.d0) alpha = 1.d0
  Eopt = energyIn + alpha*slope + 0.5d0*alpha**2*curv
  
  psi2 = psi2in + alpha*(psi2out-psi2in)
  
  H0C = H0Cin + alpha*(H0Cout-H0Cin)
  HC = 2.d0*Eopt - H0C
  
  write(*,*) '    slope ',slope
  write(*,*) '    step  ',alpha
  write(*,*) '    Eopt  ',Eopt
  
  ! check convergence
  if (abs(slope/Eopt) <= critODA) then
     converged = .true.
  else
     converged = .false.
  end if
  
  return
end subroutine oda


subroutine  convergence_test_norm(Ng,psi2in,psi2out,converged,error)
  use criteria
  
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: Ng
  real(dp), dimension(Ng), intent(in) :: psi2in,psi2out
  logical, intent(out) ::  converged
  real(dp), intent(out) :: error
  
  error = sum(abs(psi2out-psi2in))
  
  converged = .false.
  if (error <= critODA) converged=.true.
  
  return
end subroutine convergence_test_norm


!**********************************************************************
!
! Routines related to the calculation of the Hamiltonian
!
!**********************************************************************

subroutine makeH(n,Ng,psi2,H)
  !
  !  Construct non-linear Hamiltonian H(psi) 
  !
  use system1D
  use basis_set1D
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n,Ng
  real(dp), dimension(Ng), intent(in) :: psi2
  real(dp), dimension(0:n,0:n), intent(out) :: H
  
  integer :: i,ii,j

  ! Add non-linear part
  do j=0,n
     do i=j,n
        ! Linear part of the Hamiltonian
        H(i,j) = H0(i,j)
        ! Add non-linear part
        do ii=1,Ng
           H(i,j) = H(i,j) + expo(ii)*phi(i,ii)*phi(j,ii)*lambda*psi2(ii)
        end do
        H(j,i) = H(i,j)
     end do
  end do

  return
end subroutine makeH


subroutine evalH0(n,Ng,C1,C2,C1H0C2)
  !
  ! Evaluates <C1| H_0 |C2>
  !
  use system1D

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n,Ng
  real(dp), dimension(0:n), intent(in) :: C1,C2
  real(dp), intent(out) :: C1H0C2
  
  real(dp), dimension(0:n) :: H0C2

  H0C2 = matmul(H0,C2)
  
  C1H0C2 = dot_product(C1,H0C2)
  
  return
end subroutine evalH0


subroutine evalH(n,Ng,psi2,C1,C2,C1HC2)
  !
  ! Evaluates <C1| H(psi) |C2>
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: n,Ng
  real(dp), dimension(Ng), intent(in) :: psi2
  real(dp), dimension(0:n), intent(in) :: C1,C2
  real(dp), intent(out) :: C1HC2
  
  real(dp), dimension(0:n,0:n) :: H
  real(dp), dimension(0:n) :: HC2
  
  ! construct Hamiltonian matrix
  call makeH(n,Ng,psi2,H)

  ! calculate dot product
  HC2 = matmul(H,C2)

  C1HC2 = dot_product(C1,HC2)
  
  return
end subroutine evalH


subroutine initH(Ng,n)
  !
  ! Initialize Hamiltonian matrix H_0 (linear part of the Hamiltonian)
  ! and arrays phi and expo used for the calculation of the non-linear part
  !
  use system1D
  use basis_set1D
  implicit none
  integer, parameter :: dp=kind(1.d0)
  real(dp), parameter :: pi=3.1415926535897932d0

  integer, intent(in) :: Ng,n

  integer :: i,ii,j
  
  real(dp), dimension(Ng) :: zeta,w,phit,V0

  real(dp), external :: potentialV0

  phi = 0.d0

  call GaussHermite(Ng,zeta,w)

  expo = exp(zeta**2)*w*sqrt(0.5d0)

  zeta = zeta*sqrt(0.5d0)     

  do i=0,n
     ii = sf*i
     call harmonic(ii,Ng,zeta,phit)
     phi(i,:) = phit
  end do
  
  H0 = 0.d0

  ! Harmonic oscillator energy
  do i=0,n
     ii = sf*i
     H0(i,i) = real(ii,dp)+0.5d0
  end do

  ! Additional potential
  do ii=1,Ng
     V0(ii) = potentialV0(zeta(ii))
  end do

  ! Transform potential to spectral basis and add to Hamiltonian
  do j=0,n
     do i=j,n
        do ii=1,Ng
           H0(i,j) = H0(i,j) + expo(ii)*phi(i,ii)*phi(j,ii)*V0(ii)
        end do
        H0(j,i) = H0(i,j)
     end do
  end do
  
  return
end subroutine initH


!**********************************************************************
! 
! Routines to transform between basis set and grid
!
! C.M. Dion & E. Cances, Phys. Rev. E 67, 046706 (2003)
!
!**********************************************************************

subroutine basis2grid(n,Ng,c,psi)
  !
  ! Calculate wave function psi from coefficients c
  !
  use basis_set1D
  implicit none
  integer, parameter :: dp=kind(1.d0)

  integer, intent(in) :: n,Ng
  real(dp), dimension(0:n), intent(in) :: c
  real(dp), dimension(Ng), intent(out) :: psi

  integer :: i,ii

  psi = 0.d0

  do ii=1,Ng
     do i=0,n
        psi(ii) = psi(ii)+phi(i,ii)*c(i)
     end do
  end do

  return
end subroutine basis2grid