From f25e8fe3c31de6b923049be54c8deb48beaa2143 Mon Sep 17 00:00:00 2001 From: Guillaume Raffy Date: Mon, 17 Apr 2023 14:17:09 +0200 Subject: [PATCH] initial version --- cs_adiabats_2.f | 390 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 390 insertions(+) create mode 100644 cs_adiabats_2.f diff --git a/cs_adiabats_2.f b/cs_adiabats_2.f new file mode 100644 index 0000000..7b922e5 --- /dev/null +++ b/cs_adiabats_2.f @@ -0,0 +1,390 @@ + 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=50 ! 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 :: 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) + 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) + 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 + 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 +! 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 + + 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) + + 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) + + if (Etot.ge.Ebarrier(Jtot+1,m,j)) then + proba(j)=1d0 +! if (Einf(Jtot+1,m,j).eq.Elev(inist)) then + 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 + 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 + + crossJ(Jtot+1,i)=openi(Jtot+1)*crossJ(Jtot+1,i)*(2d0*Jtot+1d0) + 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(26,*)" energy jbi jai taui jbf jaf + & tauf CS" + do i=1,nlevuniq + 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 + + + + + + +************************************************************************ + +