416 lines
13 KiB
Fortran
416 lines
13 KiB
Fortran
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=119 !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))
|
|
|
|
write(6,*) 'idx_p1 = ', idx_p1
|
|
do i=np_f,1,-1 ! read in reverse to start from Rmax
|
|
write(6,*) 'i = ', i
|
|
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
|
|
|
|
|
|
|
|
|
|
|
|
|
|
************************************************************************
|
|
|
|
|