;The actual calculation of VT transition conducted in this study is based on cuSOLVER (documents can be found in https://docs.nvidia.com/cuda/cusolver/index.html).
;However, cuSOLVER is a closed-source package and requires pre-installation. It is also hardware-dependent and not entirely cross-platform.
;The program written for our calculaion is tuned speifically according to the hardware/platform and is not designed to be general.
;Here we offer a standalone IDL program to demonstrate the numerical procedures involved in the quantum calculation of transition probability (not thermal averaged)

pro trans_prob_demo
 common vv, vij,m,a0

 mydir='d:/liang/quantum/'
 
; Define constans and calcualte vibrational enengy levels
 h=6.62607d-34
 hh=h/2./!pi
 omega=fltarr(25)
 omg00=(3353.+41.2)*1.38065d-23/hh
 xe=20.6/(3353.+41.2)
 for i=0, 24 do begin
  omega[i]=omg00*(i+.5)-xe*omg00*(i+.5)*(i+.5)
 endfor

 gamma=.5
 mpm=30.0056*28.0134/(30.0056+28.0134)*1.66054d-27
 mu=7.00335*1.66054d-27
 L=1.5d-11  ; L=15 pm here
 m=gamma*gamma*mpm/mu
 a0=gamma*sqrt(hh/mu/omg00)/L

; This is the V_ij matrix calcualted accroding to equation (8), stored in a sav file 
 restore, mydir+'vij_L15.sav'

 kn=replicate(0.d,25)
 openkn=intarr(25)
 P=replicate(0.d,25)
 
 v0=5000. ; intial ion velocity in m/s

;Though in theory the transition probability is independent of the magnitude of potential, in pratical calculation it should be scaled to ensure the tuning point is within the computation regime   
 vij*=mpm*v0*v0/hh/omg00
      
 starti=5  ; Initial vibrational state
 totE0=mpm*v0*v0/hh/omg00+2.*omega[starti]/omg00
 maxkn=0
 for i=0, 24 do begin
   t1=totE0-2.*omega[i]/omg00
   if t1 gt 0. then begin
     kn[maxkn]=sqrt(m*t1)
     openkn[maxkn]=i
     maxkn++
   endif else break
 endfor
 
 if maxkn gt 1 then begin
  solveP, kn,maxkn,openkn,starti,P
  totp=total(P[0:maxkn-1])  ; sumn of all transition probabilities, should be close to 1
  print,'Sum of transition probabilities: ',totp
; printing all the open-channel states and the transtion probabilities
  for n=0, maxkn-1 do begin   
   print,FORMAT='("P(", I0, "-", I0,"): ", G)', starti,openkn[n],P[n]
  endfor
 endif

end

pro solveP, kn,maxkn,openkn,starti,P
 COMPILE_OPT idl2
 common vv, vij,m,a0

; Adjusting the x-grid range and resolution according to kn for better numerical accuracy 
 dx=(!pi/kn[0]/10.d)<0.1d
 nl=ceil(19.d/a0*kn[maxkn-1]/2.d/!pi)
 xrg=nl*2.d*!pi/kn[maxkn-1]
 
 numx=fix(xrg/2./dx)*2L+1
 x=indgen(numx)*dx-ceil(6.d/a0)

 Acoef=replicate(0.d,numx*maxkn*2,numx*maxkn*2)
 Bcoef=dblarr(numx*maxkn*2)

 kcos=dblarr(numx,maxkn)
 ksin=dblarr(numx,maxkn)
 for i=0, numx-1 do begin
   for j=0, maxkn-1 do begin
     t1=(kn[j]*i*dx) mod (2.d*!pi)
     kcos[i,j]=cos(t1)
     ksin[i,j]=sin(t1)
   endfor
 endfor
 
 for j=0, maxkn-1 do begin
 for l=0L, numx-1 do begin
  t2=exp(-a0*x[l])
  for jj=0, maxkn-1 do begin
   if vij[j,jj] eq 0. then continue
   t1=m*t2*vij[j,jj]/kn[jj]*dx

   ll=0L
   lll=l-ll
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1*17./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1*17./48.
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1*17./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1*17./48.

   ll=numx-1L
   lll=ll-l
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1*17./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1*17./48.
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1*17./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1*17./48.

   ll=1L
   lll=abs(l-ll)
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1*59./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1*59./48.
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1*59./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1*59./48.

   ll=numx-2L
   lll=abs(ll-l)
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1*59./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1*59./48.
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1*59./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1*59./48.

   ll=2L
   lll=abs(l-ll)
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1*43./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1*43./48.
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1*43./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1*43./48.

   ll=numx-3L
   lll=abs(ll-l)
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1*43./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1*43./48.
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1*43./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1*43./48.

   ll=3L
   lll=abs(l-ll)
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1*49./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1*49./48.
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1*49./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1*49./48.

   ll=numx-4L
   lll=abs(ll-l)
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1*49./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1*49./48.
   Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1*49./48.
   Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1*49./48.
   
   for ll=4L,numx-5L do begin
    lll=abs(ll-l)
    Acoef[2L*(ll+jj*numx),2L*(l+j*numx)]=ksin[lll,jj]*t1
    Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)]=kcos[lll,jj]*t1
    Acoef[2L*(ll+jj*numx),2L*(l+j*numx)+1]=-kcos[lll,jj]*t1
    Acoef[2L*(ll+jj*numx)+1,2L*(l+j*numx)+1]=ksin[lll,jj]*t1
   endfor
  endfor
  
  Acoef[2L*(l+j*numx),2L*(l+j*numx)]-=1.
  Acoef[2L*(l+j*numx)+1,2L*(l+j*numx)+1]-=1.
  
  Bcoef[2L*(l+j*numx)]=-t2*kcos[l,starti]*vij[starti,j]
  Bcoef[2L*(l+j*numx)+1]=t2*ksin[l,starti]*vij[starti,j]
 endfor
 endfor

 ;Linear-equation system solver. In this standalone IDL program we use the built-in LU-decomposion functions, it can be replaced by cuSOLVER API (once installed) or any other linear equation solver for better efficiency.
 Aludc=Acoef
 LA_LUDC, Aludc, Index,/double
 result = LA_LUSOL(Aludc, Index, Bcoef,/double)
 Ccoef = LA_LUMPROVE(TEMPORARY(Acoef), TEMPORARY(Aludc), TEMPORARY(Index), TEMPORARY(Bcoef), TEMPORARY(result),/double)
 
 Fr=dblarr(numx,maxkn)
 Fi=dblarr(numx,maxkn)
 for j=0, maxkn-1 do begin  
 for l=0, numx-1 do begin
  Fr[l,j]=Ccoef[2L*(l+j*numx)]
  Fi[l,j]=Ccoef[2L*(l+j*numx)+1]
 endfor
 endfor
 
 l=0
 Rr=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])*17./48.
 Ri=-(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])*17./48.

 l=numx-1
 Rr+=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])*17./48.
 Ri-=(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])*17./48.

 l=1
 Rr+=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])*59./48.
 Ri-=(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])*59./48.

 l=numx-2
 Rr+=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])*59./48.
 Ri-=(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])*59./48.

 l=2
 Rr+=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])*43./48.
 Ri-=(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])*43./48.

 l=numx-3
 Rr+=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])*43./48.
 Ri-=(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])*43./48.

 l=3
 Rr+=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])*49./48.
 Ri-=(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])*49./48.

 l=numx-4
 Rr+=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])*49./48.
 Ri-=(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])*49./48.
 
 for l=4, numx-5 do begin
  Rr+=(Fi[l,*]*kcos[l,*]-Fr[l,*]*ksin[l,*])
  Ri-=(Fr[l,*]*kcos[l,*]+Fi[l,*]*ksin[l,*])
 endfor

 P[0:maxkn-1]=(Rr*Rr+Ri*Ri)/kn[0:maxkn-1]/kn[starti]*m*dx*m*dx

 end