work done on 17/04/2023

.gitignore

@@ -0,0 +1,2 @@
+a.out
+cross_*.dat

cs_adiabats_2.f

@@ -17,7 +17,7 @@
       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=50        ! number of energies
+       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
@@ -25,6 +25,7 @@
 !      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
@@ -60,7 +61,10 @@
       do m=1,npar  ! loop over parities
         idx_p1=Jtot+m*1000
 
-      open(idx_p1)
+      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
@@ -71,6 +75,8 @@
 
       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)
@@ -164,7 +170,7 @@ c      read(5,*) jai, taui, jbi         ! read initial state. ja=H2O, jb=CO
       inist=0
       write(6,*) 'nlevuniq = ', nlevuniq
       do i=1,nlevuniq
-        write(6,*) 'i = ', i    
+        ! 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
@@ -181,6 +187,7 @@ c        print*, "initial state does not exist"
 !      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')
 
@@ -194,6 +201,7 @@ c        print*, "initial state does not exist"
 
       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
@@ -202,6 +210,7 @@ c        print*, "initial state does not exist"
       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
 
@@ -213,7 +222,7 @@ c        print*, "initial state does not exist"
       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
 
@@ -229,12 +238,19 @@ c        print*, "initial state does not exist"
         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
@@ -326,7 +342,10 @@ c        print*, "initial state does not exist"
 !          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)
@@ -335,9 +354,13 @@ c        print*, "initial state does not exist"
         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

readme.md

@@ -0,0 +1,16 @@
+cs_adiabats_2 provided by Jérôme Loreau
+
+count adiabats and provide cross sections and rate coefficients for asymetric top-linear molecules
+
+inputs: generated from hibridon's job.eadiab and renamed and adapted
+	`fort.<symmetry>00<jtot>` :
+		symmetry: up to 4
+		- 1 first symmetry
+		- 2 second symmetry
+		jtot: total angular momentum number
+	lvls.dat : energy levels of the systems
+		31   2   2   2   0   0   136.563
+		<index>   <qnum1>   <qnum2>   <dummy1>  <summy2>    <qnum3>    <energy>
+
+ouputs:
+	cross_<qnum1>_<qnum2>__<qnum3>.dat