crosadiab/cs_adiabats_2.f

      program findmax
! J. Loreau
! Goal: compute cross sections based on statistical method.
! proba = 0 or 1 if channel open or closed.
! Takes into account the centrifugal barrier in the adiabats.
! Developed for adiabats of COb- + COa PES
! version called from script, generates cross sections for all j_i
****** To check: np_f, file rot_levels.dat, nlevtot, mu (red mass), max number
*                adiabats

      implicit none

      integer, parameter :: np_f=102  !number of points for function, input
!      integer, parameter :: nps=10000  !number of points we work with in the splint subroutine.
      integer, parameter :: nlevtot=176  ! number of levels, some appear multiple times
      integer, parameter :: npar=2  ! number of parities from molscat
      integer, parameter :: Jmin=0  ! Jtot min
      integer, parameter :: Jmax=1  ! Jtot max
!      integer, parameter :: inist=1  ! initial state, # refers to levels.dat
       integer, parameter :: nE=5        ! number of energies
!      integer, parameter :: nf=10  !number of functions in input
!      double precision, parameter :: c = 137.035999     !in au,
!      c=1/alpha
!      double precision, parameter :: c = 299792458  !in m/s
!      double precision, parameter :: amu=1.66d-27  !in kg
!      double precision, parameter :: kb=1.38d-23 !in J/K
!      double precision, parameter :: kb2=8.617d-05 !in eV/K
      integer :: iostat
      integer :: i,j,k,m,idx,nf,nlevuniq,Jtot,idx_p1,idx_p2
      integer :: l1,l2,l1min,l1max,l2min,l2max,iE,inist,openi(Jmax+1)
      character(2) :: x1,x2,x3
      double precision :: de,pi,bid,mu,wavek2
      double precision :: R(np_f),Deltaf(np_f)
!      double precision :: func(np_f,nf)
      double precision :: Rm(5),fm(5) ! 5=max number of maxima here.
      double precision :: Rmax,fmax,ener(nE),Emin,Emax
      double precision, allocatable :: func(:,:),proba(:)
      double precision :: probalev(nlevtot,npar),probalevtot(nlevtot)
      double precision :: probaJ(Jmax+1,nlevtot)  ! Jmax+1 as index J=0 won't work
      double precision :: E_lev(nlevtot),sigidx(nlevtot),erel
      integer :: jCO(nlevtot),jH2O(nlevtot),ja(nlevtot),jb(nlevtot)
      integer :: tau(nlevtot),tau2(nlevtot),taui,taumax,tauf
      integer :: jamax,jbmax,jai,jaf,jbi,jbf,deltakr
      double precision :: Elev(nlevtot),numlev(nlevtot) ! this is too large, but to be safe
      double precision :: openadiab(npar),openadiabtot(Jmax+1)
      double precision :: crossJ(Jmax+1,nlevtot),crosstot(nlevtot)
      double precision :: Etot, Einf(Jmax+1,npar,15000)   ! 500 = max number of adiabats in file
      double precision :: Ebarrier(Jmax+1,npar,15000)   ! 500 = max number of adiabats in file
      integer :: n_adiab(Jmax+1,npar)
      character*20 name_CS
      character(len=8) :: fmt ! format descriptor
      pi=dacos(-1d0)

      fmt = '(I2.2)' ! an integer of width 2 with zeros at the left

******* read functions, extract energy of centrifugal barrier ****
      Ebarrier=0d0
      do Jtot=Jmin,Jmax       ! begin loop over Jtot
      write(400,*) "Jtot = ", Jtot

      do m=1,npar  ! loop over parities
        idx_p1=Jtot+m*1000

      open(idx_p1,status='old')

      write(6,*) 'idx_p1=', idx_p1

      read(idx_p1,*) nf
      n_adiab(Jtot+1,m)=nf    ! store the number of functions for given J, parity
      if (nf.eq.0) then  ! do nothing, only for J=0 in some cases
      else

      allocate(func(np_f,nf))
!      allocate(proba(nf))

      do i=np_f,1,-1  ! read in reverse to start from Rmax
        read(idx_p1,*) R(i),(func(i,j),j=1,nf)
!        write(6,*) 'R(', i, ') = ', R(i)
!        write(6,*) 'func( ', i,',nf=', nf,') = ', func(i,nf)
      enddo
!      rewind(idx_p1)
      close(idx_p1)
!      do i=1,nlevuniq
!       print*, numlev(i),Elev(i)
!      enddo

      do j=1,nf      ! loop over functions
      Rmax=0d0
      fmax=0d0
      Deltaf=0d0
!      idx=0
      do i=3,np_f  ! loop over distances, start from Rmax
        Deltaf(i)=func(i,j)-func(i-1,j)
!        print*, i,Deltaf(i)
        if ((Deltaf(i) .lt. 0d0) .and. (Deltaf(i-1)) .gt. 0d0) then !find maxima
!          idx=idx+1
!          Rmax=(R(i-1)+R(i))/2
!          fmax=(func(i-1,1)+func(i,1))/2
!           print*, j,R(i-1),func(i-1,j)
           if (func(i-1,j) .gt. fmax) then
             fmax=func(i-1,j)
             Rmax=R(i-1)
           endif
        endif
      enddo     ! end loop over distances
      Ebarrier(Jtot+1,m,j)=fmax            ! store max value
!      print*, Ebarrier(Jtot+1,m,j)
      Einf(Jtot+1,m,j)=func(1,j)
!      print*, Rmax,fmax
      if ((Rmax.eq.0d0).and.(Deltaf(3).lt.0d0)) then
        Ebarrier(Jtot+1,m,j)=func(1,j)       ! attractive -> barrier=E_asy
        write(400,*) "a", j,  "attractive"
      elseif ((Rmax.eq.0d0).and.(Deltaf(3).gt.0d0)) then
        Ebarrier(Jtot+1,m,j)=1d10         ! arbitrary large value, proba=0
        write(400,*) j,  "repulsive"
      elseif (fmax.lt.func(1,j)) then   
        Ebarrier(Jtot+1,m,j)=func(1,j)    ! takes care of local min in the well
        write(400,*) "b", j,  "attractive"
      endif
* note: the order in the "if" is important! Do not invert 2 and 3.

      write(400,*) "Ebarrier = ", Ebarrier(Jtot+1,m,j)

      enddo     ! end loop over functions
      deallocate(func)

      endif  ! end if over nf=0 (to treat the J=0 case)

      enddo   ! end loop over parities m

      enddo     !enf loop over Jtot

!      print*, Ebarrier(1,1,:)

*******
** we now know Ebarrier and Einf for each J, parity, # adiabat
** we can start the loop over the energies.


c      read(5,*) jai, taui, jbi         ! read initial state. ja=H2O, jb=CO

      open(9,file='lvls.dat',status='old')
      read(9,*)         ! determine the number of unique levels, nlevuniq
      nlevuniq=0
      do i=1,nlevtot
        write (6,*) 'i = ', i
        read(9,*) sigidx(i),jH2O(i),tau2(i),bid, bid,jCO(i),E_lev(i)  
        write (6,*) 'sigidx = ', sigidx(i)          
        if (sigidx(i).ne.sigidx(i-1)) then
          nlevuniq=nlevuniq+1
          numlev(nlevuniq)=sigidx(i)
          Elev(nlevuniq)=E_lev(i)
          ja(nlevuniq)=jH2O(i)
          jb(nlevuniq)=jCO(i)
          tau(nlevuniq)=tau2(i)
!          print*,nlevuniq,tau(nlevuniq)
        endif
      enddo
      close(9)
      jamax=maxval(ja)
      jbmax=maxval(jb)

      do jai=0,9
      taumax=2*jai+1
!      print*, "taumax=",taumax
      do taui=1,taumax
      do jbi=0,7
!      print*, jai,taui,jbi

      inist=0
      write(6,*) 'nlevuniq = ', nlevuniq
      do i=1,nlevuniq
        ! write(6,*) 'i = ', i    
        if ((ja(i).eq.jai) .and. (jb(i).eq.jbi) .and. (tau(i).eq.taui))
     &  then
        inist=i    ! find number of level corresponding to initial state
        endif
      enddo
      print*, "inist=", inist, "with energy", Elev(inist)
      if (inist .ne. 0) then  ! if inist =0, don't do anything (wrong parity)
c        print*, "initial state does not exist"
!      print*, "inist=", inist

      write (x1,fmt) ja(inist) ! converting integer to string using a 'internal file'
      write (x2,fmt) tau(inist) !
      write (x3,fmt) jb(inist) !
!      print*,x1,x2,x3
      name_CS='cross_'//x1//'_'//x2//'__'//x3//'.dat'
!      print*, name_CS
      write(6, *) 'opening file ', name_CS
      open(26,file=name_CS,status='replace')
      open(400,file='output_adiab',status='replace')






      Emin=1d0
      Emax=1d3

      do iE=1,nE    ! loop over energies
      ener(iE)=Emin+(iE-1)*(Emax-Emin)/(nE-1)  !define scale for energy
       write(6,*) 'ener(', iE, ') = ', ener(iE)
!      ener(1)=500d0
      write(400,*) 'energy: ',ener(iE)
!      ener(iE)=10+(iE-1)*20d0  !define scale for energy
!      ener(iE)=ener(iE)-Elev(inist)     ! energy in cm-1

      probaJ=0d0
      crosstot=0d0

      write(6, *) 'Jmin=', Jmin, ' Jmax=', Jmax
      do Jtot=Jmin,Jmax       ! begin loop over Jtot
!      write(6,*) "E = ", ener(iE), "Jtot = ", Jtot

      openi(Jtot+1)=0
      probalev=0d0
      probalevtot=0d0
      openadiab=0d0

      do m=1,npar  ! loop over parities

      nf=n_adiab(Jtot+1,m)
      write(6,*) 'nf = ', nf
      if (nf.eq.0) then  ! do nothing, only for J=0 for some cases
      else

      allocate(proba(nf))

      do j=1,nf         ! loop over adiabats

        do i=1,nlevuniq ! determine the level associated with the adiabat
!          if (Einf(Jtot+1,m,j).eq.Elev(i)) then
          if (abs(Einf(Jtot+1,m,j)-Elev(i)).lt.1d-4) then
            k=i
          endif
        enddo

        Etot=ener(iE)+Elev(inist)
        write(6, *) 'Elev(', inist, ') = ', Elev(inist)
        write(6, *) 'iE = ', iE

        write(6, *) 'Etot = ', Etot,
     &' Ebarrier(Jtot+1,m,j)=', Ebarrier(Jtot+1,m,j)
        if (Etot.ge.Ebarrier(Jtot+1,m,j)) then
          proba(j)=1d0
!          if (Einf(Jtot+1,m,j).eq.Elev(inist)) then
          write(6, *) 'Einf = ', Einf(Jtot+1,m,j),
     & ' Elev =', Elev(inist)
          if (abs(Einf(Jtot+1,m,j)-Elev(inist)).lt.1d-4) then
            openi(Jtot+1)=openi(Jtot+1)+1   ! count the number of open initial channels 
            stop
          endif
        else
          proba(j)=0d0
        endif


        probalev(k,m)=probalev(k,m)+proba(j)  ! sum the probas for each unique level and parity
        openadiab(m)=openadiab(m)+proba(j)   ! counts the number of open channels per parity
      enddo     ! end loop over adiabats
      write(400,*) "Jtot = ", Jtot, "open channels for Jtot=", Jtot, 
     &        "and parity: ", m, "= ", openadiab(m)

!      print*, "J = ", Jtot
!      do i=1,nlevuniq
!        print*, i,probalev(i,m),probalev(i,m)/sum(proba)
!      enddo

      deallocate(proba)
!      Ptot=sum(proba)
!      print*, Ptot
!      do i=1,nlevuniq ! determine the level associated with the adiabat
!        if (func(1,j).eq.Elev(i)) then
!        k=i
!      endif
!      enddo
!
!      do j=1,nf
!        if k(j)=i then
!        probalev=probalev+proba(j)
!        endif

      endif  ! end if over nf=0 (to treat the J=0 case)

      enddo   ! end loop over parities m
      openadiabtot(Jtot+1)=sum(openadiab)       ! sum over m
      write(400,*) "channel n°,       open channels,      probability:"
      do i=1,nlevuniq
      do m=1,npar
        probalevtot(i)=probalevtot(i)+probalev(i,m)
      enddo
!        print*, i,probalev(i,m),probalev(i,m)/sum(proba)

      if (openadiabtot(Jtot+1).eq.0d0) then
        probaJ(Jtot+1,i)=0d0    ! to avoid infinity below
        write(400,*) i,probalevtot(i),probaJ(Jtot+1,i)
      else
        probaJ(Jtot+1,i)=probalevtot(i)/sum(openadiab)
        write(400,*) i,probalevtot(i),probaJ(Jtot+1,i)
      endif
      enddo


      enddo     !end loop over Jtot


!      call cross(probaJ,)

**** Compute the cross section based on S-matrix (= probas)

      mu=10.956d0 ! reduced mass
      wavek2=2d0*mu*1822.8885d0*ener(iE)/219474.63d0  ! in a.u.
      crossJ=0d0
!      jbi=jb(inist)  ! already defined
!      jai=0  ! COa
      do Jtot=Jmin,Jmax
!         print*, Jtot, openi(Jtot+1)
!        print*, openadiabtot(Jtot+1)
        if (openadiabtot(Jtot+1).eq.0d0) then
          crossJ(Jtot+1,:)=0d0
        else
        do i=1,nlevuniq
          jaf=ja(i)
          jbf=jb(i)
          tauf=tau(i)

* this sum is not necessary, already taken into account by adiabats #
!          l1min=abs(Jtot-ji)       ! change  ji and jf here!
!          l1max=Jtot+ji
!          l2min=abs(Jtot-jf)
!          l2max=Jtot+jf
!          print*, 'l1...', l1min, l1max, l2min, l2max
!          do l1=l1min,l1max
!          do l2=l2min,l2max
            crossJ(Jtot+1,i)=crossJ(Jtot+1,i)+
!     &      (deltakr(ji,jf)*deltakr(l1,l2)-probaJ(Jtot+1,i))**2
!     &      (deltakr(ji,jf)-probaJ(Jtot+1,i))**2  !no square!
     &      abs(deltakr(jbi,jbf)*deltakr(jai,jaf)*deltakr(taui,tauf)
     &          -probaJ(Jtot+1,i))
!          enddo
!          enddo

          write(6, *) "openi(", Jtot+1, ")=", openi(Jtot+1)

          crossJ(Jtot+1,i)=openi(Jtot+1)*crossJ(Jtot+1,i)*(2d0*Jtot+1d0)
          write(6, *) "crossJ(", Jtot+1, ", ", i,")=", crossJ(Jtot+1,i)
          crosstot(i)=crosstot(i)+crossJ(Jtot+1,i)
          if (crossJ(Jtot+1,i).ne.0d0) then
            write(300+i,*) ener(iE), Jtot, crossJ(Jtot+1,i)
          endif
        enddo
        endif
      enddo


      write(6, *) "writing header"
      write(26,*)"    energy       jbi   jai   taui    jbf   jaf 
     & tauf      CS"
      write(6, *) 'nlevuniq = ', nlevuniq
      do i=1,nlevuniq
        write(6, *) "crosstot(", i, ")=", crosstot(i)
        crosstot(i)=crosstot(i)*pi/((2d0*jai+1d0)*(2d0*jbi+1d0)*wavek2)
        crosstot(i)=crosstot(i)*((5.291772d-1)**2)
        if (crosstot(i).ne.0d0) then
        write(26,111) ener(iE),jbi,jai,taui,jb(i),ja(i),tau(i),
     &                crosstot(i)
        endif
      enddo

      enddo   ! end loop over energies
      close(400)
      close(26)

c      endif ! end if over inist=0
c      enddo  ! end loops over jai,jbi,tau
c      enddo
c      enddo


!      print*, probaJ
111   format(1(ES14.6),6I6,ES14.4)
112   format(17(E16.8,1h,))



      endif ! end if over inist=0
      enddo  ! end loops over jai,jbi,tau
      enddo
      enddo

       end program 


       function deltakr(a,b)
       implicit none
       integer :: a,b,deltakr
       if (a.eq.b) then
         deltakr=1d0
       else
         deltakr=0d0
       endif
       end function






************************************************************************