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=142 !number of points for function, input ! integer, parameter :: nps=10000 !number of points we work with in the splint subroutine. integer, parameter :: nlevtot=165 ! 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=5 ! Jtot max ! integer, parameter :: inist=1 ! initial state, # refers to levels.dat integer, parameter :: nE=500 ! 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 wont 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='level.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,5 taumax=2*jai+1 ! print*, "taumax=",taumax do taui=1,taumax do jbi=0,10 ! 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, dont 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=8d2 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-3) 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 ************************************************************************