gpoda1dg.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
!
! Grid method in 1D
!
! Claude M. Dion and Eric Cances
!
! Version 1.0 (February 2007)
! 
!**********************************************************************

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

module system1Dg
  real(kind(1.d0)) :: mass,lambda
  real(kind(1.d0)), dimension(:), allocatable :: kinetic,V
end module system1Dg

module grid_desc
  real(kind(1.d0)) :: d
  real(kind(1.d0)), dimension(2) :: grid_ends
  real(kind(1.d0)), dimension(:), allocatable :: grid
end module grid_desc


program GPODA1Dg
  !
  ! main program
  !
  use system1Dg
  use grid_desc

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

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

  integer :: itMax,iter
  real(dp) :: mu,error,H0psiin,Hinpsiin
  
  logical :: guess_from_file,converged

  real(dp), external :: integrate1D

  write(*,*) "GPODA1Dg"
  write(*,*)

  !
  ! Initialization
  !
  
  call read_params(ng,itMax,guess_from_file)
 
  write(*,*) "Initialization"

  allocate(grid(ng),kinetic(ng),V(ng))
  call init_linear(ng)
  call fft_init(ng)

  allocate(psi2in(ng),psi2out(ng),psi(ng))

  write(*,*)
  if(guess_from_file) then 
     ! read initial guess from file guess1Dg.data
     write(*,*) "  Read initial guess from file guess1Dg.data"
     call read_psi(ng,"guess1Dg.data",psi)
     psi2in = psi**2
  else
     ! use ground state of H_0 as guess 
     write(*,*) "  Compute the ground state of H_0"
     call guess_gauss(ng,psi) ! start from Gaussian
     psi2in = 0.d0
     call linear_gs(ng,psi2in,0.0D0,psi,mu,psi,psi2in)
  endif
  write(*,*)

  ! compute initial values of H_0 psi and H(psi) psi
  call evalH0(ng,psi,psi,H0psiin)
  call evalH(ng,psi2in,psi,psi,Hinpsiin)

  !
  ! Iterations of the ODA
  !
  do iter=1,ItMax
     write(*,*)
     write(*,*) "Iteration ",iter
     write(*,*) 
     write(*,*) "  Compute the ground state of H(psi_in)"
     call linear_gs(ng,psi2in,0.D0,psi,mu,psi,psi2out)
     write(*,*) "  Optimal damping"
     call oda(ng,psi,psi2in,psi2out,H0psiin,Hinpsiin,psi2in, &
          & H0psiin,Hinpsiin,converged)
     write(*,*) 
     if(converged) exit
  end do

  if(converged) then
     write(*,*) 
     write(*,*) 'Convergence achieved in',iter,' iterations'
     write(*,*) ' --> mu =',mu
     write(*,*) ' --> E = ',0.5d0*(H0psiin+Hinpsiin)
     write(*,*)
     
     write(*,*)
     write(*,*) "Checking self-consistency"
     write(*,*)         
     call linear_gs(ng,psi2in,0.0D0,psi,mu,psi,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

  ! write out order parameter on a grid
  call write_grid(ng,psi,"gs1Dg.data")

  deallocate(grid,kinetic,V,psi2in,psi2out,psi)
  
end program GPODA1Dg


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

subroutine read_params(ng,itMax,guess_from_file)
  !  
  ! Reads parameters from the file params1Dg.in
  !
  use criteria
  use system1Dg
  use grid_desc

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

  integer, intent(out) :: ng,itMax
  logical, intent(out) :: guess_from_file

  real(dp) :: zmin,zmax
  
  namelist /params1Dg/ &
       & mass, &         ! mass of the boson (in a.u.)   
       & lambda, &       ! non-linearity (in a.u.)
       & ng, &           ! number of grid points 
       & zmin, &         ! first point in z
       & zmax, &         ! last point in z
       & critODA, &      ! convergence criterion for the ODA
       & critIP, &       ! convergence criterion for the inverse power
       & critCG, &       ! convergence criterion for the conjugated gradient
       & itMax, &        ! maximum number of iterations of the ODA
       & guess_from_file ! read initial guess from file guess1Dg.data?

  ! read parameters from file
  open(unit=8,file="params1Dg.in",status='old')
  read(8,params1Dg)
  close(8)
  
  ! check validity of some parameters
  if(ng < 1) then
     write(*,*) "Invalid input parameter ng =",ng
     stop
  end if

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

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

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

  ! transfer coordinate-dependent values to an array
  grid_ends(1) = zmin
  grid_ends(2) = zmax

  ! print out parameters 
  write(*,90) mass,lambda,ng,grid_ends

90 format("Parameters:"/,/," mass =      ",e16.8,/,&
        & " non-linearity = ",e16.8,/,/, &
        & " Number of grid points: ",i4,/ &
        & " Extent in z: ",e14.6," -- ",e14.6)

  return
end subroutine read_params


subroutine read_psi(ng,infile,psi)
  !
  ! Read initial order parameter
  !
  ! *** It is assumed that the grid of the input file is the same
  ! *** as that defined in the parameter file
  !
  use system1Dg
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: ng
  character(len=*), intent(in) :: infile
  real(dp), dimension(ng), intent(out) :: psi

  integer :: k
  real(dp) :: z

  ! Read order parameter
  open(13,file=infile,status='old')

  do k = 1,ng
     read(13,*) z,psi(k)
  end do
  
66 close(13)

  ! Renormalize psi
  call normalize(ng,psi)

  return
end subroutine read_psi


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

  integer :: k

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

  do k = 1,ng
     write(12,88) grid(k),psi(k)
88   format(e14.6,1x,e16.8)
  end do

  close(12)

  return
end subroutine write_grid


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

subroutine linear_gs(ng,psi2,shift,guess,eigenval,psiout,psi2out)
  !
  ! Calculate the ground state of H(psi) 
  ! by the inverse power method
  !
  use criteria
  use grid_desc

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: ng
  real(dp), dimension(ng), intent(in) :: psi2,guess
  real(dp), intent(in) :: shift
  real(dp), intent(out) :: eigenval
  real(dp), dimension(ng), intent(out) :: psiout,psi2out
  
  integer :: iter
  real(dp) :: eigval_old,eigval_new
  real(dp), dimension(ng) :: u,v

  real(dp), external :: integrate1D

  ! guess for the ground state  
  u = guess
  
  ! first guess for the eigenvalue should be below ground state value
  eigval_old = -1.d20
  
  ! initial guess for v
  v = 0.d0

  ! inverse power loop
  iter = 0  
  do
     iter = iter+1
     call solveH(ng,psi2,u,v)
     eigval_new = 1.d0/integrate1D(ng,d,u*v)
     u = v/sqrt(integrate1D(ng,d,v**2))
     ! convergence achieved?
     if (abs(eigval_new-eigval_old) <= critIP) exit
     ! no, iterate
     eigval_old = eigval_new
     v = u/eigval_old
  end do
  
  ! convergence achieved
  eigenval = eigval_new+shift
  write(*,88) iter
88 format("   Inverse Power converged in ",i6," iterations")
  write(*,*) "   --> mu =",eigenval

  psiout = u
  psi2out = psiout**2

  return
end subroutine linear_gs


subroutine solveH(ng,psi2,u,v)
  !
  ! Solve H(psi).v = u for v using a conjugated gradient method
  !
  ! On input, v contains a guess for its value
  !
  use criteria
  use grid_desc

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: ng
  real(dp), dimension(ng), intent(in) :: psi2,u
  real(dp), dimension(ng), intent(inout) :: v
  
  real(dp) :: alpha,beta,num,den
  real(dp), dimension(ng) :: Ad,r,dd,Qr

  real(dp), external :: integrate1D

  call operatorH(ng,psi2,v,Ad)
  r = u-Ad
  dd = r
  Qr = r
  num = integrate1D(ng,d,Qr*r) 

  do
     if (sqrt(integrate1D(ng,d,r**2)) <= critCG) exit
     call operatorH(ng,psi2,dd,Ad)
     den = integrate1D(ng,d,Ad*dd)
     alpha = num/den
     v = v + alpha*dd
     r = r - alpha*Ad
     Qr = r
     den = num
     num = integrate1D(ng,d,Qr*r) 
     beta = num/den
     dd = r + beta*dd
  end do

  return
end subroutine solveH


subroutine oda(ng,psiout,psi2in,psi2out,H0psiin,Hinpsiin,psi2, &
     & H0psi,Hpsi,converged)
  !
  ! Optimal Damping Algorithm
  !
  use criteria
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: ng
  real(dp), dimension(ng), intent(in) :: psiout,psi2in,psi2out
  real(dp), intent(in) :: H0psiin,Hinpsiin
  real(dp), dimension(ng), intent(out) :: psi2
  real(dp), intent(out) :: H0psi,Hpsi
  logical, intent(out) :: converged

  real(dp) :: H0psiout,Houtpsiout,Hinpsiout
  real(dp) :: energyIn,slope,curv,alpha,Eopt
  
  call evalH0(ng,psiout,psiout,H0psiout)
  call evalH(ng,psi2out,psiout,psiout,Houtpsiout)
  call evalH(ng,psi2in,psiout,psiout,Hinpsiout)
  
  energyIn = 0.5d0*(H0psiin+Hinpsiin)
  slope = Hinpsiout - Hinpsiin
  curv = Hinpsiin + Houtpsiout - 2.d0*Hinpsiout+H0psiout-H0psiin
  
  alpha = -slope/curv
  if(alpha >= 1.d0) alpha = 1.d0
  Eopt = energyIn + alpha*slope + 0.5d0*alpha**2*curv
  
  psi2 = psi2in + alpha*(psi2out-psi2in)
  
  H0psi = H0psiin + alpha*(H0psiout-H0psiin)
  Hpsi = 2.d0*Eopt - H0psi
  
  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
  use grid_desc
  
  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 operatorH(ng,psi2,psiin,Hpsiin)
  !
  ! Applies the non-linear Hamiltonian H(psi) to the order parameter psiin
  ! (resulting order parameter is returned in Hpsiin)
  !
  use system1Dg

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: ng
  real(dp), dimension(ng), intent(in) :: psi2,psiin
  real(dp), dimension(ng), intent(out) :: Hpsiin

  real(dp), dimension(ng) :: psitmp

  ! Calculate kinetic energy

  ! transform to momentum space
  call fourier1D(ng,psiin,psitmp,+1)

  ! kinetic energy
  psitmp = kinetic*psitmp

  ! transform back to position space
  call fourier1D(ng,psitmp,Hpsiin,-1)

  ! Calculate potential energy and non-linearity
  Hpsiin = Hpsiin + (V+lambda*psi2)*psiin

  return
end subroutine operatorH


subroutine evalH0(ng,psiA,psiB,psiAH0psiB)
  !
  ! Evaluates <psiA| H_0 |psiB>
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: ng
  real(dp), dimension(ng), intent(in) :: psiA,psiB
  real(dp), intent(out) :: psiAH0psiB
  
  real(dp), dimension(ng) ::  psi2
  
  ! Set psi to zero and evaluate <psiA| H(psi2) |psiB>
  psi2 = 0.d0
  call evalH(ng,psi2,psiA,psiB,psiAH0psiB)
  
  return
end subroutine evalH0


subroutine evalH(ng,psi2,psiA,psiB,psiAHpsiB)
  !
  ! Evaluates <psiA| H(psi) |psiB>
  !
  use grid_desc
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: ng
  real(dp), dimension(ng), intent(in) :: psi2,psiA,psiB
  real(dp), intent(out) :: psiAHpsiB
  
  real(dp), dimension(ng) :: HpsiB
  real(dp), external :: integrate1D
  
  call operatorH(ng,psi2,psiB,HpsiB)
  psiAHpsiB = integrate1D(ng,d,psiA*HpsiB)
  
  return
end subroutine evalH


subroutine init_linear(ng)
  !
  ! Initialize arrays for the linear part of the Hamiltonian
  !
  use system1Dg
  use grid_desc
  implicit none
  integer, parameter :: dp=kind(1.d0)
  real(dp), parameter :: pi=3.1415926535897932d0

  integer, intent(in) :: ng

  integer :: i,k,imax
  real(dp) :: fact

  real(dp), external :: potentialV
  
  ! construct grid
  d = (grid_ends(2)-grid_ends(1))/real(ng-1,dp)
  do i = 1,ng
     grid(i) = grid_ends(1) + real(i-1,dp)*d
  end do

  ! calculate kinetic operator in momentum space
  fact = 2.d0*(pi/(real(ng,dp)*d))**2/mass
  kinetic(1) = 0.d0

  if(2*(ng/2) == ng) then
     imax = ng/2
     kinetic(ng) = fact*real(ng/2,dp)**2
  else
     imax = (ng+1)/2
  end if

  do i = 2,imax
     kinetic(2*i-2) = fact*real(i-1,dp)**2
     kinetic(2*i-1) = kinetic(2*i-2)
  end do

  ! calculate potential
  do k = 1,ng
     V(k) = potentialV(grid(k))
  end do

  return
end subroutine init_linear


!**********************************************************************
!
! Utility routines
!
!**********************************************************************

subroutine guess_gauss(ng,psi)
  !
  ! Construct Gaussian as initial guess
  !
  use grid_desc
  implicit none
  integer, parameter :: dp=kind(1.d0)
  real(dp), parameter :: epsilon = 1d-8
  
  integer, intent(in) :: ng
  real(dp), dimension(ng), intent(out) :: psi

  integer :: k
  real(dp) :: centre,sigma2

  ! center of grid
  centre = (grid_ends(2)-grid_ends(1))*0.5d0 + grid_ends(1)

  ! width
  sigma2 = (grid_ends(2) - centre)**2/(-2.d0*log(epsilon))

  ! Gaussian
  do k = 1,ng
     psi(k) = exp(-(grid(k)-centre)**2/(2.d0*sigma2))
  end do

  call normalize(ng,psi)

  return
end subroutine guess_gauss


subroutine normalize(ng,psi)
  !
  ! Renormalizes order parameter psi
  !
  use grid_desc
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, intent(in) :: ng
  real(dp), dimension(ng), intent(inout) :: psi

  real(dp) :: norm
  real(dp), dimension(ng) :: psi2
  real(dp), external :: integrate1D
  
  psi2 = psi**2
  norm = integrate1D(ng,d,psi2)
  psi = psi/sqrt(norm)

  return
end subroutine normalize