gpoda3dg.f90

Ground state of the time-independent Gross-Pitaevskii equation

     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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(in) :: psi2,u
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(inout) :: v
  
  real(dp) :: alpha,beta,num,den
  real(dp), dimension(ng(1),ng(2),ng(3)) :: Ad,r,dd,Qr

  real(dp), external :: integrate3D
  
  call operatorH(ng,psi2,v,Ad)
  r = u-Ad
  dd = r
  Qr = r
  num = integrate3D(ng(1),ng(2),ng(3),d(1),d(2),d(3),Qr*r) 

  do
     if(sqrt(integrate3D(ng(1),ng(2),ng(3),d(1),d(2),d(3),r**2)) <= critCG) exit
     call operatorH(ng,psi2,dd,Ad)
     den = integrate3D(ng(1),ng(2),ng(3),d(1),d(2),d(3),Ad*dd)
     alpha = num/den
     v = v + alpha*dd
     r = r - alpha*Ad
     Qr = r
     den = num
     num = integrate3D(ng(1),ng(2),ng(3),d(1),d(2),d(3),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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(in) :: psiout,psi2in,psi2out
  real(dp), intent(in) :: H0psiin,Hinpsiin
  real(dp), dimension(ng(1),ng(2),ng(3)), 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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), 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 system3Dg

  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,psiin
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(out) :: Hpsiin

  integer :: i,j,k
  real(dp) :: ktmp
  real(dp), dimension(ng(1),ng(2),ng(3)) :: psitmp

  ! Calculate kinetic energy
  
  ! transform to momentum space
  call fourier3D(ng,psiin,psitmp,+1)

  ! kinetic energy
  do k = 1,ng(3)
     do j = 1,ng(2)
        ktmp = kinetic(j,2)+kinetic(k,3)
        do i = 1,ng(1)
           psitmp(i,j,k) = (kinetic(i,1)+ktmp)*psitmp(i,j,k)
        end do
     end do
  end do

  ! transform back to position space
  call fourier3D(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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(in) :: psiA,psiB
  real(dp), intent(out) :: psiAH0psiB
  
  real(dp), dimension(ng(1),ng(2),ng(3)) ::  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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(in) :: psi2,psiA,psiB
  real(dp), intent(out) :: psiAHpsiB
  
  real(dp), dimension(ng(1),ng(2),ng(3)) :: HpsiB
  real(dp), external :: integrate3D
  
  Call operatorH(ng,psi2,psiB,HpsiB)
  psiAHpsiB = integrate3D(ng(1),ng(2),ng(3),d(1),d(2),d(3),psiA*HpsiB)
  
  return
end subroutine evalH


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

  integer, dimension(0:3), intent(in) :: ng

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

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

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

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

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

  ! calculate potential
  do k = 1,ng(3)
     do j = 1,ng(2)
        do i = 1,ng(1)
           V(i,j,k) = potentialV(grid(i,1),grid(j,2),grid(k,3))
        end do
     end do
  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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(out) :: psi

  integer :: i,j,k
  real(dp) :: tmp,tmpz
  real(dp), dimension(3) :: centre,sigma2

  ! center of grid
  do i = 1,3
     centre(i) = (grid_ends(2,i)-grid_ends(1,i))*0.5d0 + grid_ends(1,i)
  end do

  ! widths
  do i = 1,3
     sigma2(i) = (grid_ends(2,i) - centre(i))**2/(-2.d0*log(epsilon))
  end do

  ! Gaussian
  do k = 1,ng(3)  
     tmpz = exp(-(grid(k,3)-centre(3))**2/(2.d0*sigma2(3)))
     do j = 1,ng(2)
        tmp = exp(-(grid(j,2)-centre(2))**2/(2.d0*sigma2(2))) * tmpz
        do i = 1,ng(1)
           psi(i,j,k) = exp(-(grid(i,1)-centre(1))**2/(2.d0*sigma2(1))) * tmp
        end do
     end do
  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, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(inout) :: psi

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

  return
end subroutine normalize