From 46471a8c12c7021b4236952d1265e6415100e81d Mon Sep 17 00:00:00 2001 From: Guillaume Raffy Date: Mon, 5 Jun 2023 14:19:23 +0200 Subject: [PATCH] work done on 17/04/2023 --- .gitignore | 2 ++ cs_adiabats_2.f | 31 +++++++++++++++++++++++++++---- readme.md | 16 ++++++++++++++++ 3 files changed, 45 insertions(+), 4 deletions(-) create mode 100644 .gitignore create mode 100644 readme.md diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4b2fa8e --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +a.out +cross_*.dat diff --git a/cs_adiabats_2.f b/cs_adiabats_2.f index 7b922e5..865690e 100644 --- a/cs_adiabats_2.f +++ b/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 diff --git a/readme.md b/readme.md new file mode 100644 index 0000000..7bf7d11 --- /dev/null +++ b/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.00` : + 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 + + +ouputs: + cross____.dat \ No newline at end of file