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

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

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

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


program GPODA3Ds
  !
  ! main program
  !
  use basis_set3D
  use system3D

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3) :: 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(*,*) "GPODA3Ds"
  write(*,*)

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

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

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

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

  write(*,*)
  if(guess_from_file) then 
     ! read initial guess from file guess3Ds.data
     write(*,*) "  Read initial guess from file guess3Ds.data"
     call read_coeff(n,"guess3Ds.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,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)*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
     
  ! convention: lowest coefficient is positive
  if (C(0,0,0) < 0.d0) C = -C
  
  ! write out coefficients of the ground state
  call write_coeff(n,C,"gs3Ds.data")
  
  ! write out order parameter on a grid
  if(output_grid) call write_grid(n,C,"gs3Ds_grid.data")

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


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

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

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

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

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

  ! read parameters from file
  open(unit=8,file="params3Ds.in",status='old')
  read(8,params3Ds)
  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_y < 1) then
     write(*,*) "Invalid input parameter n_y =",n_y
     stop
  elseif(n_y > 91) then
     write(*,*) "Too many basis functions n_y =",n_y
     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_y
  n(3) = n_z
  n(0) = maxval(n(1:3))
  symmetric(1) = symmetric_x
  symmetric(2) = symmetric_y
  symmetric(3) = symmetric_z

  ! calculate constants to reconstruct HO number from index
  do k=1,3
     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:3) = 2*sf*n(1:3)+1
  Ng(0) = maxval(Ng(1:3))

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

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

  if(any(symmetric)) then
     sym = "Symmetric in"
     if(symmetric(1)) sym = trim(sym)//" x"
     if(symmetric(2)) sym = trim(sym)//" y"
     if(symmetric(3)) 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 system3D
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: n
  character(len=*), intent(in) :: infile
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(out) :: c

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

     do l=1,3
        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),i(3)) /= 0.d0) then
        write(*,*) 'Warning!  Component ',i,' defined more than once,', &
             & ' last value kept'
     end if

     c(i(1),i(2),i(3)) = 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 system3D
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: n
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(in) :: c
  character(len=*), intent(in) :: outfile

  integer :: i,j,k,ii,jj,kk

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

  do i=0,n(1)
     ii = sf(1)*i
     do j=0,n(2)
        jj = sf(2)*j
        do k=0,n(3)
           kk = sf(3)*k
           write(12,*) ii,jj,kk,C(i,j,k)
        enddo
     enddo
  enddo

  close(12)

  return
end subroutine write_coeff


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

  integer :: ng_x,ng_y,ng_z,i,j,k,ii,jj,kk
  real(dp) :: xmin,xmax,ymin,ymax,zmin,zmax,dx,dy,dz
  real(dp), dimension(:), allocatable :: x,y,z,phit
  real(dp), dimension(:,:), allocatable :: phi
  real(dp), dimension(:,:,:), allocatable :: psi,ccz,cyz

  namelist /grid3D/ &
       & 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

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

  ! construct grid
  allocate(x(ng_x),y(ng_y),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

  dy = (ymax-ymin)/real(ng_y-1,dp)
  do jj = 1,ng_y
     y(jj) = ymin + real(jj-1,dp)*dy
  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(3),ng_z),phit(ng_z))
  phi = 0.d0  
  do k = 0,n(3)
     kk = sf(3)*k
     call harmonic(kk,ng_z,z,phit)
     phi(k,:) = phit
  end do
  deallocate(phit)  

  allocate(ccz(0:n(1),0:n(2),ng_z))
  ccz = 0.d0
  do kk = 1,ng_z
     do k = 0,n(3)
        ccz(:,:,kk) = ccz(:,:,kk)+phi(k,kk)*C(:,:,k)
     end do
  end do
  deallocate(phi)
  
  ! transform y coordinate
  allocate(phi(0:n(2),ng_y),phit(ng_y))
  phi = 0.d0

  do j = 0,n(2)
     jj = sf(2)*j
     call harmonic(jj,ng_y,y,phit)
     phi(j,:) = phit
  end do
  deallocate(phit)

  allocate(cyz(0:n(1),ng_y,ng_z))
  cyz = 0.d0
  do jj = 1,ng_y
     do j = 0,n(2)
        cyz(:,jj,:) = cyz(:,jj,:)+phi(j,jj)*ccz(:,j,:)
     end do
  end do 
  deallocate(phi,ccz)

  ! 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_y,ng_z))
  psi = 0.d0
  do ii = 1,ng_x
     do i = 0,n(1)
        psi(ii,:,:) = psi(ii,:,:)+phi(i,ii)*cyz(i,:,:)
     end do
  end do
  deallocate(phi,cyz)