SUBROUTINE basis_functions ! !========================================================================================= ! Written By: J. Crawford ! ! This routine obtains the asymptotic diatomic eigenfunctions as a function of ! Jacobi mass-scaled s. It also maps the functions to the APH grid and stores ! the functions in a binary file to be used later. ! ! Variables: ! jac_function = array storing Jacobi basis functions ! j_in = diatomic angular momentum j in initial wave packet ! jmax = max value of diatomic angular momentum j ! jtotal = value of total angular momentum ! l_in = space fixed orbital angular momentum of initial wave packet ! nvib_in = diatomic vibrational state in initial wave packet ! parity = parity of initial wave packet ! req_chanls = number of unique arrangement channels ! sdim = dimension of small s Jacobi grid (read in from namelist) ! s_dim = dimension of small s Jacobi grid (read in from binary) ! s_max = max small s Jacobi grid value ! s_min = min small s JAcobi grid value ! s_vals = array of small s Jacobi grid values !========================================================================================= Use Numeric_Kinds_Module USE CommonInfo_Module USE SurfaceJacobi_Module USE SurfaceDelves_Module USE SurfaceAph_Module USE QuantumNumber_Module USE Print_Module USE SymGroup_Module USE APH_Module IMPLICIT NONE !======================================================================================== ! I N P U T !========================================================================================= ! I N T E R N A L S CHARACTER(LEN=1) :: c1 CHARACTER(LEN=2) :: j2 CHARACTER(LEN=100):: filename INTEGER :: i, ichanl, ivib, ijrot, istat, jcheck, ccheck, neig, sdim, s_dim, tcdim REAL(dp) :: s_min, s_max !========================================================================================= ! A L L O C A T A B L E REAL(dp),ALLOCATABLE :: jac_function(:,:), s_vals(:) REAL(dp),ALLOCATABLE :: wfrot(:), wfvib(:) !COMPLEX(dp),ALLOCATABLE :: wfrot(:), wfvib(:) !========================================================================================= ! F U N C T I O N S CHARACTER(LEN=1) :: label1 CHARACTER(LEN=2) :: label2 INTEGER :: aphindex, aphindex_tc save !========================================================================================= ! N A M E L I S T S NAMELIST / surface_grid / sdim, td_dim !========================================================================================= ! Read in NAMELIST parameters OPEN(UNIT=Input_Unit,FILE=Trim(inputfile1),IOSTAT=istat,STATUS='old',ACTION='read') IF (istat.ne.0) STOP 'OPEN failed - surface_function' READ(Input_Unit,surface_grid) CLOSE(UNIT=Input_Unit,IOSTAT=istat,STATUS='keep') IF (istat.ne.0) STOP 'CLOSE failed - surface_function' !========================================================================================= ! Check input parameter validity IF (j_in.gt.jmax) jmax=j_in !========================================================================================= ! Output pertinent info to screen WRITE(*,21) 'jtotal ',jtotal WRITE(*,21) 'jmax ',jmax WRITE(*,21) 'nvib_in',nvib_in WRITE(*,21) 'j_in ',j_in WRITE(*,21) 'l_in ',l_in WRITE(*,21) 'parity ',parity WRITE(*,22) 'Calculating Jacobi Basis Functions' 21 FORMAT(1X,a,8('.'),I3) 22 FORMAT(1X,a) !========================================================================================= ! Determine number of unique arrangement channels nchanls=3 req_chanls=nchanls IF (system.eq.'AAA') req_chanls=1 IF (system.eq.'ABB') req_chanls=Min(2,nchanls) !========================================================================================== ! Calculate Jacobi rovib functions tcdim=ntheta*nchi firstcallj = .true. firstcalld = .true. !========================================================================================= ! Obtain rovibrational eigenfunctions in Jacobi and Delves coordinates ALLOCATE( wfrot(tcdim), wfvib(tcdim) ) DO ichanl=1,nchanls ! DO (1) : channel loop c1=label1(ichanl) CALL jacobi_grid(ichanl) DO ijrot = 0, jmax ! DO (2) : j loop j2=label2(ijrot) IF (ichanl.le.req_chanls) THEN ! IF (1) : check for req_chanls IF (.not.use_pvs.or.ijrot.eq.j_in) THEN ! IF (2) : check for use_pvs or j_in !================================================================================== ! Obtain vibrational basis functions for j=ijrot in arrangement channel=ichanl CALL jacobi_basis( ijrot, ichanl, sdim ) !CALL delves_surface( ijrot, ichanl) ENDIF ! IF (2) : use_pvs or j_in ENDIF ! IF (1) : req_chanls !=================================================================================== ! Read in rovibrational small s rovib functions calculated in jacobi_surface routine IF (ichanl.gt.req_chanls) c1=label1(req_chanls) filename='jacobi_surface_c'//c1//'_'//j2//'.bin' OPEN(UNIT=bin_unit1,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown') READ(bin_unit1) ccheck READ(bin_unit1) jcheck IF (jcheck.ne.ijrot) STOP 'wrong jjval readin : surface_functions' READ(bin_unit1) neig READ(bin_unit1) s_dim !IF (sdim_read.ne.sdim) STOP 'wrong sdim readin : surface_functions' 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 IF(ALLOCATED(jac_function)) DEALLOCATE(jac_function) ALLOCATE( jac_function(s_dim,neig)) READ(bin_unit1) jac_function ! IF (ichanl.eq.nchanls) THEN ! IF (ichanl.ne.1.and.ijrot.ne.j_in) THEN ! CLOSE(UNIT=bin_unit1, STATUS='delete') ! ELSE ! CLOSE(bin_unit1) ! ENDIF ! ELSE CLOSE(bin_unit1) ! ENDIF IF (.not.use_pvs) THEN ! IF (3) : check for use_pvs !=========================================================================================== ! Open output file to store Jacobi small s cap theta rovib functions until needed in analysis c1=label1(ichanl) filename='jacobi_totsurface_c'//c1//'_'//j2//'.bin' OPEN(UNIT=bin_unit0,FILE=TRIM(BinOutdir)//TRIM(filename),FORM='unformatted',STATUS='unknown') !============================================================================ ! Obtain rotational component of asymptotic basis function mapped to APH grid IF (ntheta.eq.1) THEN wfrot=0.d0 !wfrot=Cmplx(1.d0,0.d0,dp) ELSE CALL rotational(ijrot,0,tcdim, wfrot,'O','SF') ! Obtain cap theta legendre polynomial ENDIF DO ivib=1,neig ! DO (3) : vibrational loop !====================================================== ! Obtain vibrational basis function mapped to APH grid CALL vibrational(s_dim, s_vals, s_min, s_max, jac_function(:,ivib),tcdim,wfvib,'O') wfvib=wfvib*wfrot WRITE(bin_unit0) wfvib ENDDO ! DO (3) : vib loop CLOSE(bin_unit0) ENDIF ! IF (3) use_pvs ENDDO ! DO (2) : j loop firstcallj=.true. ENDDO ! DO (1) : channel loop !========================================================================================= ! Deallocate arrays DEALLOCATE(wfrot,wfvib) END SUBROUTINE basis_functions