gpoda3dg.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 3D
!
! Claude M. Dion and Eric Cances
!
! Version 1.0 (February 2007)
! 
!**********************************************************************

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

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

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


program GPODA3Dg
  !
  ! main program
  !
  use system3Dg
  use grid_desc

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

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

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

  write(*,*) "GPODA3Dg"
  write(*,*)

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

  allocate(grid(ng(0),3),kinetic(ng(0),3),V(ng(1),ng(2),ng(3)))
  call init_linear(ng)
  call fft_init(ng)

  allocate(psi2in(ng(1),ng(2),ng(3)),psi2out(ng(1),ng(2),ng(3)), &
       & psi(ng(1),ng(2),ng(3)))

  write(*,*)
  if(guess_from_file) then 
     ! read initial guess from file guess3Dg.data
     write(*,*) "  Read initial guess from file guess3Dg.data"
     call read_psi(ng,"guess3Dg.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(1)*ng(2)*ng(3),error/real(ng(1)*ng(2)*ng(3),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,"gs3Dg.data")

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


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

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

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

  integer, dimension(0:3), intent(out) :: ng
  integer, intent(out) :: itMax
  logical, intent(out) :: guess_from_file

  integer :: ng_x,ng_y,ng_z
  real(dp) :: xmin,xmax,ymin,ymax,zmin,zmax
  
  namelist /params3Dg/ &
       & mass, &         ! mass of the boson (in a.u.)   
       & lambda, &       ! non-linearity (in a.u.)
       & ng_x, &         ! number of grid points in x
       & ng_y, &         ! number of grid points in y
       & ng_z, &         ! number of grid points in z
       & xmin, &         ! first point in x
       & xmax, &         ! last point in x
       & ymin, &         ! first point in y
       & ymax, &         ! last point in y
       & 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 guess3Dg.data?

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

  if(ng_y < 1) then
     write(*,*) "Invalid input parameter ng_y =",ng_y
     stop
  end if

  if(ng_z < 1) then
     write(*,*) "Invalid input parameter ng_z =",ng_z
     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 arrays
  ng(1) = ng_x
  ng(2) = ng_y
  ng(3) = ng_z
  ng(0) = maxval(ng(1:3))
  grid_ends(1,1) = xmin
  grid_ends(2,1) = xmax
  grid_ends(1,2) = ymin
  grid_ends(2,2) = ymax
  grid_ends(1,3) = zmin
  grid_ends(2,3) = zmax

  ! print out parameters 
  write(*,90) mass,lambda,ng(1:3),ng(1)*ng(2)*ng(3),grid_ends

90 format("Parameters:"/,/," mass =      ",e16.8,/,&
        & " non-linearity = ",e16.8,/,/, &
        & " Number of grid points: ",i4," x ",i4," x ",i4," = ",i10,/ &
        & " Extent in x: ",e14.6," -- ",e14.6,/ &
        & "           y: ",e14.6," -- ",e14.6,/ &
        & "           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 system3Dg
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: ng
  character(len=*), intent(in) :: infile
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(out) :: psi

  integer :: i,j,k
  real(dp) :: x,y,z

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

  do k = 1,ng(3)  
     do j = 1,ng(2)
        do i = 1,ng(1)
           read(13,*) x,y,z,psi(i,j,k)
        end do
     end do
  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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(in) :: psi
  character(len=*), intent(in) :: outfile

  integer :: i,j,k

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

  do k = 1,ng(3)  
     do j = 1,ng(2)
        do i = 1,ng(1)
           write(12,88) grid(i,1),grid(j,2),grid(k,3),psi(i,j,k)
88         format(3(e14.6,1x),e16.8)
        end do
     end do
  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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(in) :: psi2,guess
  real(dp), intent(in) :: shift
  real(dp), intent(out) :: eigenval
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(out) :: psiout,psi2out
  
  integer :: iter
  real(dp) :: eigval_old,eigval_new
  real(dp), dimension(ng(1),ng(2),ng(3)) :: u,v

  real(dp), external :: integrate3D

  ! 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/integrate3D(ng(1),ng(2),ng(3),d(1),d(2),d(3),u*v)
     u = v/sqrt(integrate3D(ng(1),ng(2),ng(3),d(1),d(2),d(3),v**2))
     ! convergence achieved?
     if (abs(eigval_new-eigval_old) <= critIP) exit
     ! no, iterate
     eigval_old = eigval_new