gpoda2ds.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 2D
!
! Claude M. Dion and Eric Cances
!
! Version 1.0 (January 2007)
! 
!**********************************************************************

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

module system2D
  integer, dimension(2) :: sf
  real(kind(1.d0)) :: wxwz,lambda
  real(kind(1.d0)), dimension(:,:), allocatable :: V0,ener
  logical, dimension(2) :: symmetric
end module system2D

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


program GPODA2Ds
  !
  ! main program
  !
  use basis_set2D
  use system2D

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

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

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

  write(*,*) "GPODA2Ds"
  write(*,*)

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

  allocate(ener(0:n(1),0:n(2)))
  call init_ener(n)

  allocate(phi(0:n(0),Ng(0),2),expo(Ng(1),Ng(2)),V0(Ng(1),Ng(2)))
  call init_nonHO(Ng,n)

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

  write(*,*)
  if(guess_from_file) then 
     ! read initial guess from file guess2Ds.data
     write(*,*) "  Read initial guess from file guess2Ds.data"
     call read_coeff(n,"guess2Ds.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,0) = 1.d0
     psi2in = 0.d0
     call linear_gs(n,Ng,psi2in,0.0D0,C,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 linear_gs(n,Ng,psi2in,0.D0,C,mu,C,psi2out)
     write(*,*) "  Optimal damping"
     call oda(n,Ng,C,psi2in,psi2out,H0Cin,HinCin,psi2in,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 linear_gs(n,Ng,psi2in,0.0D0,C,mu,C,psi2out)
     call convergence_test_norm(Ng,psi2in,psi2out,converged,error)
     write(*,89) error,Ng(1)*Ng(2),error/real(Ng(1)*Ng(2),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) < 0.d0) C = -C
  
  ! write out coefficients of the ground state
  call write_coeff(n,C,"gs2Ds.data")

  ! write out order parameter on a grid
  if(output_grid) call write_grid(n,C,"gs2Ds_grid.data")  

  deallocate(ener,phi,expo,V0,psi2in,psi2out,C)
  
end program GPODA2Ds


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

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

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

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

  integer :: k,n_x,n_z
  logical :: symmetric_x,symmetric_z
  character(len=16) :: sym
  
  namelist /params2Ds/ &
       & lambda, &          ! non-linear factor
       & wxwz, &            ! trap frequency ratio omega_x/omega_z
       & n_x, &             ! highest basis function in x
       & n_z, &             ! highest basis function in z
       & symmetric_x, &     ! symmetric in x?
       & symmetric_z, &     ! symmetric 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 guess2Ds.data?
       & output_grid        ! ouput order parameter on a grid?

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

  if(n_z < 1) then
     write(*,*) "Invalid input parameter n_z =",n_z
     stop
  elseif(n_z > 91) then
     write(*,*) "Too many basis functions n_z =",n_z
     write(*,*) "Maximum value is 91"
     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
  n(1) = n_x
  n(2) = n_z
  n(0) = maxval(n(1:2))
  symmetric(1) = symmetric_x
  symmetric(2) = symmetric_z

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

  ! calculate number of spatial grid points
  Ng(1:2) = 2*sf*n(1:2)+1
  Ng(0) = maxval(Ng(1:2))

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

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

  if(any(symmetric)) then
     sym = "Symmetric in"
     if(symmetric(1)) sym = trim(sym)//" x"
     if(symmetric(2)) sym = trim(sym)//" z"
     write(*,*) sym
  end if
        
  write(*,*)

  return
end subroutine read_params


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

  integer :: l
  integer, dimension(2) :: 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(1),i(2),tmp

     do l=1,2
        if(symmetric(l)) then
           if(2*(i(l)/2) /= i(l)) then
              write(*,*) 'Component ',i,' of input file ignored'
              cycle reading
           end if
           i(l) = i(l)/sf(l)
        end if
        if(i(l) > n(l) .or. i(l) < 0) then
           write(*,*) 'Component ',i,' of input file ignored'
           cycle reading
        end if
     end do

     if(c(i(1),i(2)) /= 0.d0) then
        write(*,*) 'Warning!  Component ',i,' defined more than once,', &
             & ' last value kept'
     end if

     c(i(1),i(2)) = tmp
  end do reading

66 close(13)

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

  return
end subroutine read_coeff


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

  integer :: i,k,ii,kk

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

  do i=0,n(1)
     ii = sf(1)*i
     do k=0,n(2)
        kk = sf(2)*k
        write(12,*) ii,kk,C(i,k)
     end do
  end do

  close(12)

  return
end subroutine write_coeff


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

  integer :: ng_x,ng_z,i,k,ii,kk
  real(dp) :: xmin,xmax,zmin,zmax,dx,dz
  real(dp), dimension(:), allocatable :: x,z,phit
  real(dp), dimension(:,:), allocatable :: phi
  real(dp), dimension(:,:), allocatable :: psi,cz

  namelist /grid2D/ &
       & ng_x, & ! number of grid points in x
       & ng_z, & ! number of grid points in z
       & xmin, & ! first point in x
       & xmax, & ! last point in x
       & zmin, & ! first point in z
       & zmax    ! last point in z

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

  ! construct grid
  allocate(x(ng_x),z(ng_z))

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

  dz = (zmax-zmin)/real(ng_z-1,dp)
  do kk = 1,ng_z
     z(kk) = zmin + real(kk-1,dp)*dz
  end do


  ! transform z coordinate
  allocate(phi(0:n(2),ng_z),phit(ng_z))
  phi = 0.d0  
  do k = 0,n(2)
     kk = sf(2)*k
     call harmonic(kk,ng_z,z,phit)
     phi(k,:) = phit
  end do
  deallocate(phit)  

  allocate(cz(0:n(1),ng_z))
  cz = 0.d0
  do kk = 1,ng_z
     do k = 0,n(2)
        cz(:,kk) = cz(:,kk)+phi(k,kk)*C(:,k)
     end do
  end do
  deallocate(phi)
  
  ! transform x coordinate
  allocate(phi(0:n(1),ng_x),phit(ng_x))
  phi = 0.d0

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

  ! write order parameter
  open(unit=12,file=trim(outfile))

  do kk = 1,ng_z
     do ii = 1,ng_x
        write(12,88) x(ii),z(kk),psi(ii,kk)
88      format(3(e14.6,1x),e16.8)