gpoda3ds.f90

Ground state of the time-independent Gross-Pitaevskii equation

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

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

  close(12)

  deallocate(x,y,z,psi)

  return
end subroutine write_grid


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

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

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: n,Ng
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(in) :: psi2
  real(dp), intent(in) :: shift
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(in) :: guess
  real(dp), intent(out) :: eigenval
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(out) :: C
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(out) :: psi2out
  
  integer :: iter
  real(dp) :: eigval_old,eigval_new
  real(dp), dimension(0:n(1),0:n(2),0:n(3)) :: u,v
  real(dp), dimension(Ng(1),Ng(2),Ng(3)) :: psi

  ! 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(n,Ng,psi2,u,v)
     eigval_new = 1.d0/sum(u*v)
     u = v/sqrt(sum(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

  C = u

  call basis2grid(n,Ng,C,psi)
  psi2out = psi**2

  return
end subroutine linear_gs


subroutine solveH(n,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

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: n,Ng
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(in) :: psi2
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(in) :: u
  real(dp), dimension(0:n(1),0:n(2),0:n(3)),intent(inout) :: v
  
  real(dp) :: alpha,beta,num,den
  real(dp), dimension(0:n(1),0:n(2),0:n(3)) :: Ad,r,d,Qr
  
  call operatorH(n,Ng,psi2,v,Ad)
  r = u-Ad
  d = r
  Qr = r
  num = sum(Qr*r)

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

  return
end subroutine solveH


subroutine oda(n,Ng,Cout,psi2in,psi2out,H0Cin,HinCin,psi2,H0C,HC,converged)
  !
  ! Optimal Damping Algorithm
  !
  use criteria
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: n,Ng
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(in) :: Cout
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(in) :: psi2in,psi2out
  real(dp), intent(in) :: H0Cin,HinCin
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(out) :: psi2
  real(dp), intent(out) :: H0C,HC
  logical, intent(out) :: converged

  real(dp) :: H0Cout,HoutCout,HinCout
  real(dp) :: energyIn,slope,curv,alpha,Eopt
  
  call evalH0(n,Ng,Cout,Cout,H0Cout)
  call evalH(n,Ng,psi2out,Cout,Cout,HoutCout)
  call evalH(n,Ng,psi2in,Cout,Cout,HinCout)
  
  energyIn = 0.5d0*(H0Cin+HinCin)
  slope = HinCout - HinCin
  curv = HinCin + HoutCout - 2.d0*HinCout+H0Cout-H0Cin
  
  alpha = -slope/curv
  if(alpha >= 1.d0) alpha = 1.d0
  Eopt = energyIn + alpha*slope + 0.5d0*alpha**2*curv
  
  psi2 = psi2in + alpha*(psi2out-psi2in)
  
  H0C = H0Cin + alpha*(H0Cout-H0Cin)
  HC = 2.d0*Eopt - H0C
  
  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
  
  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(n,Ng,psi2,C,HC)
  !
  ! Applies the non-linear Hamiltonian H(psi) to the state vector C
  ! (resulting vector is returned in HC)
  !
  use system3D

  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: n,Ng
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(in) :: psi2
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(in) :: C
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(out) :: HC
  
  real(dp), dimension(Ng(1),Ng(2),Ng(3)) :: HCx

  ! Transform from basis set to grid points

  call basis2grid(n,Ng,C,HCx)

  ! Calculate extra potential and non-linearity
  HCx = HCx*(V0+lambda*psi2)

  ! Transform back to basis set
  call grid2basis(n,Ng,HCx,HC)

  ! Add H.O. energy
  HC = HC + ener*C

  return
end subroutine operatorH


subroutine evalH0(n,Ng,C1,C2,C1H0C2)
  !
  ! Evaluates <C1| H_0 |C2>
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: n,Ng
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(in) :: C1,C2
  real(dp), intent(out) :: C1H0C2
  
  real(dp), dimension(Ng(1),Ng(2),Ng(3)) ::  psi2
  
  ! Set psi to zero and evaluate <C1| H(psi2) |C2>
  psi2 = 0.d0
  call evalH(n,Ng,psi2,C1,C2,C1H0C2)
  
  return
end subroutine evalH0


subroutine evalH(n,Ng,psi2,C1,C2,C1HC2)
  !
  ! Evaluates <C1| H(psi) |C2>
  !
  implicit none
  integer, parameter :: dp=kind(1.d0)
  
  integer, dimension(0:3), intent(in) :: n,Ng
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(in) :: psi2
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(in) :: C1,C2
  real(dp), intent(out) :: C1HC2
  
  real(dp), dimension(0:n(1),0:n(2),0:n(3)) :: HC2
  
  Call operatorH(n,Ng,psi2,C2,HC2)
  C1HC2 = sum(C1*HC2)
  
  return
end subroutine evalH


subroutine init_ener(n)
  !
  ! Initialize array ener used to calculate the energy 
  ! of the harmonic oscillator
  !
  use system3D
  implicit none
  integer, parameter :: dp=kind(1.d0)

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

  integer :: i,j,k,ii,jj,kk
  
  do k=0,n(3)
     kk = sf(3)*k
     do j=0,n(2)
        jj = sf(2)*j
        do i=0,n(1)
           ii = sf(1)*i
           ener(i,j,k) = wxwz*(real(ii,dp)+0.5d0)+wywz*(real(jj,dp)+0.5d0) &
                & +(real(kk,dp)+0.5d0)
        end do
     end do
  end do

  return
end subroutine init_ener


subroutine init_nonHO(Ng,n)
  !
  ! Initialize arrays phi, expo, and V0 used for
  ! the calculation of the non-linear part of the Hamiltonian
  !
  use system3D
  use basis_set3D
  implicit none
  integer, parameter :: dp=kind(1.d0)
  real(dp), parameter :: pi=3.1415926535897932d0

  integer, dimension(0:3), intent(in) :: Ng,n

  integer :: i,ii,jj,kk,m
  real(dp), dimension(Ng(0),3) :: zeta
  real(dp), dimension(Ng(0)) :: w,phit
  real(dp), dimension(Ng(0),3) :: etmp

  real(dp), external :: potentialV0

  phi = 0.d0

  do m=1,3
     call GaussHermite(Ng(m),zeta(1,m),w)

     etmp(:Ng(m),m) = exp(zeta(:Ng(m),m)**2)*w(:Ng(m))*sqrt(0.5d0)

     zeta(:Ng(m),m) = zeta(:Ng(m),m)*sqrt(0.5d0)     

     do i=0,n(m)
        ii = sf(m)*i
        call harmonic(ii,Ng(m),zeta(1,m),phit)
        phi(i,:Ng(m),m) = phit(:Ng(m))
     end do
  end do
  
  do kk=1,Ng(3)
     do jj=1,Ng(2)
        do ii=1,Ng(1)
           V0(ii,jj,kk) = potentialV0(zeta(ii,1),zeta(jj,2),zeta(kk,3))
        end do
     end do
  end do

  do kk=1,Ng(3)
     do jj=1,Ng(2)
        do ii=1,Ng(1)
           expo(ii,jj,kk) = etmp(ii,1)*etmp(jj,2)*etmp(kk,3)
        end do
     end do
  end do

  return
end subroutine init_nonHO


!**********************************************************************
! 
! Routines to transform between basis set and grid
!
! C.M. Dion & E. Cances, Phys. Rev. E 67, 046706 (2003)
!
!**********************************************************************

subroutine basis2grid(n,Ng,c,psi)
  !
  ! Calculate wave function psi from coefficients c
  !
  use basis_set3D
  implicit none
  integer, parameter :: dp=kind(1.d0)

  integer, dimension(0:3), intent(in) :: n,Ng
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(in) :: c
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(out) :: psi

  integer :: i,j,k,ii,jj,kk
  real(dp), dimension(0:n(1),0:n(2),Ng(3)) :: ccx
  real(dp), dimension(0:n(1),Ng(2),Ng(3)) :: cxx

  ! transform z coordinate
  ccx = 0.d0
  
  do kk=1,Ng(3)
     do k=0,n(3)
        ccx(:,:,kk) = ccx(:,:,kk)+phi(k,kk,3)*C(:,:,k)
     end do
  end do
  
  ! transform y coordinate
  cxx = 0.d0
  
  do jj=1,Ng(2)
     do j=0,n(2)
        cxx(:,jj,:) = cxx(:,jj,:)+phi(j,jj,2)*ccx(:,j,:)
     end do
  end do 
  
  ! transform x coordinate
  psi = 0.d0

  do ii=1,Ng(1)
     do i=0,n(1)
        psi(ii,:,:) = psi(ii,:,:)+phi(i,ii,1)*cxx(i,:,:)
     end do
  end do

  return
end subroutine basis2grid


subroutine grid2basis(n,Ng,psi,c)
  !
  ! Calculate coefficients c from wave function psi 
  !
  use basis_set3D
  implicit none
  integer, parameter :: dp=kind(1.d0)

  integer, dimension(0:3), intent(in) :: n,Ng
  real(dp), dimension(Ng(1),Ng(2),Ng(3)), intent(in) :: psi
  real(dp), dimension(0:n(1),0:n(2),0:n(3)), intent(out) :: c

  integer :: i,j,k,ii,jj,kk
  real(dp), dimension(0:n(1),0:n(2),Ng(3)) :: ccx
  real(dp), dimension(0:n(1),Ng(2),Ng(3)) :: cxx

  ! transform x coordinate

  cxx = 0.d0

  do i=0,n(1)
     do ii=1,Ng(1)
        cxx(i,:,:) = cxx(i,:,:)+phi(i,ii,1)*psi(ii,:,:)*expo(ii,:,:)
     end do
  end do

  ! transform y coordinate

  ccx = 0.d0

  do j=0,n(2)
     do jj=1,Ng(2)
        ccx(:,j,:) = ccx(:,j,:)+phi(j,jj,2)*cxx(:,jj,:)
     end do
  end do

  ! transform z coordinate

  c = 0.d0

  do k=0,n(3)
     do kk=1,Ng(3)
        c(:,:,k) = c(:,:,k)+phi(k,kk,3)*ccx(:,:,kk)
     end do
  end do

  return
end subroutine grid2basis