Subroutine initial_wavepacket

Use Numeric_Kinds_Module
Use QuantumNumber_Module
Use CommonInfo_Module
Use APH_Module
Use Jacobi_Module
Use Gaussian_Module
Use SurfaceJacobi_Module
USe SurfaceAph_Module
Use SymGroup_Module

Implicit None

!=========================================================================================
! This routine outputs the symmetrized initial atom-diatom wave packet in APH coordinates. 
!   The wave packet is calculated using Jacobi mass-scaled wave functions and mapped to 
!   APH coordinates and symmetrized into its component irreducible representaions.
!
! Variables:
!   phi          = array containing the symmetrized initial wave packet
!   symmat       = array containing symmetrization matrix
!   to_propagate = logical array containing information about the irreducible 
!                  representaions present in the initial wave packet.
! 
!=========================================================================================
!   I N P U T

!=========================================================================================
!   O U T P U T

!=========================================================================================
!   P A R A M E T E R S

LOGICAL,PARAMETER :: debug=.true.

!=========================================================================================
!   I N T E R N A L S

Character(len=2)   :: j2, label2
Character(len=100) :: filename
Integer     :: aphindex, ispt, irho, itheta, ichi, jchi, istat, nchir
Integer     :: whichir, indexchi, ccheck, jcheck, neig, s_dim
Real(dp)    :: jacobicoef, aphcoef, probamp, x, y, tol
Real(dp)    :: rho_sum, theta_sum, chi_sum, norm, s_min, s_max

!=========================================================================================
!   A L L O C A T A B L E
Real(dp),Allocatable    :: symmat(:,:), s_vals(:),wftemp(:)
Complex(dp),Allocatable :: phi(:), resultvec(:), splitvec(:), wftemp2(:)

!=========================================================================================


tol=1.0d-12

!=========================================================================================
! Allocate
  ALLOCATE(wftemp(ndim),phi(ndim))
  ALLOCATE(to_propagate(nirrep))

  phi=Cmplx(0.d0,0.d0,dp)
  !wftemp=Cmplx(0.d0,0.d0,dp)

!========================================================================================
!  Switch to initial Jacobi coordinates

  Call jacobi_grid(1)

!========================================================================================
!  Read in appropriate information
  j2=label2(j_in)
  filename='jacobi_surface_c1_'//j2//'.bin'
  OPEN(UNIT=bin_unit1,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
  READ(bin_unit1) ccheck
  IF (ccheck.ne.1) STOP 'wrong channel : initital_wavepacket'
  READ(bin_unit1) jcheck
  IF (jcheck.ne.j_in) STOP  'wrong j value : initial_wavepacket'
  READ(bin_unit1) neig
  READ(bin_unit1) s_dim
  IF(ALLOCATED(s_vals)) DEALLOCATE(s_vals)
  ALLOCATE( s_vals(s_dim))
  READ(bin_unit1) s_vals
  READ(bin_unit1) s_min
  READ(bin_unit1) s_max
  !CLOSE(UNIT=bin_unit1,STATUS='delete') 
  CLOSE(bin_unit1)
!========================================================================================
!  Calculate initial diatomic rotational function

  IF (ntheta.eq.1) THEN

    phi=Cmplx(1.d0,0.d0,dp)

  ELSE
    CALL rotational(j_in, l_in, ndim, wftemp, 'I', 'SF')

    phi=Cmplx(wftemp,0.d0,dp)

    wftemp=0.d0
  ENDIF

!=========================================================================================
!  Interpolate initial Diatomic Vibrational Function

  Call vibrational(s_dim, s_vals, s_min, s_max, vib_eig_in, ndim, wftemp, 'I')

  print*,wftemp(1000:1010)

  DEALLOCATE(s_vals)

  phi=phi*wftemp

  DEALLOCATE(wftemp)
  ALLOCATE(wftemp2(ndim))
!=========================================================================================
!    Calculate gaussian wavepacket 
! Keep the same for now j(initial)=0  = l(initial)

  Call S_wavepacket(wftemp2)!spacket)

  phi = phi*wftemp2

!==================================================
!    Obtain initial energy distribution (eta)

  Call eta

!==================================================
!    Combine components to obtain total wave packet

IF (debug) THEN
  filename="iwf/total_iwf.txt"
  OPEN(UNIT=output_unit0,FILE=TRIM(outdir)//TRIM(filename),IOSTAT=istat,STATUS='replace',ACTION='write')
  IF (istat.ne.0) STOP 'OPEN failed - initial_wavepacket'

  
  filename="iwf/xy_grid.txt"
  OPEN(UNIT=output_unit2,FILE=TRIM(outdir)//TRIM(filename),IOSTAT=istat,STATUS='replace',ACTION='write')
  IF (istat.ne.0) STOP 'OPEN failed - initial_wavepacket'

ENDIF

  Do irho=1,nrho
    Do itheta=1,ntheta
      Do ichi=1,nchi
         ispt=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)

         IF (ntheta.eq.1) THEN
           jacobicoef=1.d0
           aphcoef=Sqrt(rho_val(irho))
           norm=1.d0
         ELSE
           jacobicoef=1.d0/(larges(ispt)*smalls(ispt))
           aphcoef=Sqrt(rho_val(irho)**5)/4.d0
           norm=1.d0/Sqrt(2.d0) ! because aph coords cover config space twice
         ENDIF

         phi(ispt)=norm*jacobicoef*aphcoef*phi(ispt)*Sqrt(weights(itheta))
     
         IF (debug) THEN   
           IF (irho.eq.ind_rho_infty) THEN
             probamp=Real(phi(ispt))**2+Aimag(phi(ispt))**2
             IF (probamp.lt.1.d-15) probamp=0.d0
             x=TAN(0.5d0*theta_val(itheta))*COS(chi_val(ichi))
             y=TAN(0.5d0*theta_val(itheta))*SIN(chi_val(ichi))
             Write(output_unit0,*) real(x,sp),real(y,sp),real(probamp,sp)
             Write(output_unit2,*) real(x,sp),real(y,sp)
           ENDIF
         ENDIF

      End Do
    End Do
  End Do


IF (debug) THEN
  CLOSE(UNIT=output_unit0,IOSTAT=istat,STATUS='keep')
  IF (istat.ne.0) STOP 'CLOSE failed - initial_wavefunction'

  CLOSE(UNIT=output_unit2,IOSTAT=istat,STATUS='keep')
  IF (istat.ne.0) STOP 'CLOSE failed - initial_wavefunction'
ENDIF

  DEALLOCATE(wftemp2)

!=========================================================================================
!    Symmetrize

  ALLOCATE( symmat(nchi,nchi), resultvec(nchi), splitvec(nchi) )

  filename='symmat.bin'
  OPEN(111,file=TRIM(BinOutdir)//TRIM(filename),form='unformatted')
    READ(111) symmat
  CLOSE(111)

IF (debug) THEN
  filename="iwf/total_iwf_sym.txt"
  OPEN(UNIT=output_unit0,FILE=TRIM(outdir)//TRIM(filename),IOSTAT=istat,STATUS='replace',ACTION='write')
  IF (istat.ne.0) STOP 'OPEN failed - initial_wavepacket'
ENDIF

  DO irho=1,nrho
  DO itheta=1,ntheta
    DO ichi=1,nchi
      ispt=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)
      splitvec(ichi)=phi(ispt)
    ENDDO
    resultvec=Matmul(Transpose(symmat),splitvec)
    DO ichi=1,nchi
      ispt=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)
      phi(ispt)=resultvec(ichi)
      IF (debug) THEN
        IF (irho.eq.ind_rho_infty) THEN
        x=TAN(0.5d0*theta_val(itheta))*COS(chi_val(ichi))
        y=TAN(0.5d0*theta_val(itheta))*SIN(chi_val(ichi))
        probamp=Real(phi(ispt))**2+Aimag(phi(ispt))**2
        IF (probamp.lt.1.d-15) probamp=0.d0
          Write(output_unit0,*) real(x,sp),real(y,sp),real(probamp,sp)
        ENDIF
      ENDIF
    ENDDO
  ENDDO
 ENDDO

IF (debug) THEN
  CLOSE(UNIT=output_unit0,IOSTAT=istat,STATUS='keep')
  IF (istat.ne.0) STOP 'CLOSE failed - initial_wavefunction'
ENDIF

CALL normalize(ndim,nchi,phi)

!=========================================================================================
! Check which irred. reps to be propagated
to_propagate=.false.
ircount: DO whichir=1,nirrep
  nchir=nchi_ir(whichir)
  indexchi=indexchi_array(whichir)
  DO irho=1,nrho
    DO itheta=1,ntheta
      DO ichi=1,nchir
        ispt=aphindex(ntheta,nchi,irho,itheta,ichi+indexchi,.false.)
        IF (Abs(Real(phi(ispt),dp)).gt.tol) THEN
          to_propagate(whichir)=.true.
          CYCLE ircount
        ENDIF    
      ENDDO
    ENDDO
  ENDDO
ENDDO ircount 

print*,'to_propagate =',to_propagate

!=========================================================================================
! Output symmetrized wave function to binary

 filename = 'iwf_sym.bin'
 OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown')
 WRITE(bin_unit0) phi
 CLOSE(bin_unit0,STATUS='keep')

 Deallocate(phi)

End Subroutine initial_wavepacket