!Copolymer code works; March 25,2020 (started writing on Feb 29, 2020)! !------------- specify #Maximum Quanta--------------------! Module common_variables implicit none Character :: Answer Double Precision :: dw, w, dt, wavenumber, wavelength, AbmaxX, AbmaxY, DeltaEps, weV Double Precision :: Abmaxtot, Fmax Double Precision, parameter :: t=4.5 !time! Double Precision, External :: FC Double Precision, EXTERNAL :: Repulsion Double Precision, EXTERNAL :: Dopant Double Precision, EXTERNAL :: EnergyCT Double Precision, EXTERNAL :: Tefunction Double Precision, EXTERNAL :: Thfunction Double Precision, EXTERNAL :: Jfunction Double Precision, EXTERNAL :: Dipoles !Double Precision, External :: disorder_table Double Precision :: lambP, lambN, lambEP, lambEN, NPROB, rand1, rand2, rand1N, rand2N Double Precision :: randR, theta, rand1Nnew, Meandist, difference integer :: i1, j1, i2, j2, f1, f2, i11, j11, i22, j22, i, j, vp1, vn1, vp2, vn2, neh, f11, f22 integer :: i3, i33, j3, j33, v3, v33, x3, kount5, kount4, kount3, vp11, vn11, vp22, vn22 integer :: x4, v4, x5 integer :: INFO, LWORK, kount1, kount2, x1, x2, t1, v1, v2, v11, v22, Kbase, configuration integer :: whole Double precision :: wholeEl, wholemean, wholedif, standD, wholePmax Double precision :: DelJTintra, DelJTinter integer, parameter :: Realization = 1000 !realization! integer, parameter :: N = 10 integer, parameter :: Cell = 1 integer, parameter :: vibmax = 0 integer, parameter :: vibmaxT = 0 integer, parameter :: B1 = N*(vibmax+1) !1PE states)! integer, parameter :: B2 = (N*(N-1)*vibmax*(vibmax+1))/2 !2PE states! integer, parameter :: B3 = (N)*(vibmax+1)*(vibmax+2) !CTn.n states! integer, parameter :: B4 = ((N)*(N-2)*vibmaxT*(vibmaxT+1)*(vibmaxT+2))/3 !CTn.nvib states! integer, parameter :: B5 = ((N-2)*(2*N-3)*vibmax*(vibmax-1)*(2*vibmax+2))/12 !3PE states! integer, parameter :: B = (B1+B2+B3+B4+B5) integer, parameter :: LDA = B integer, parameter :: Z= 10000 integer, parameter :: y= 10 Double Precision :: beta = 2.35d0 Double Precision :: sigma = 0.5 Double Precision :: randmean = 0.0 Double Precision :: distance = 2.1d0 Double Precision :: danion = 4.1d0 Double Precision :: lamb = 0.0d0 Double Precision :: Sfactor = 1.0 Double Precision, parameter :: VMU = 2.88 Double Precision :: wdstar = 0.0 Double Precision :: wastar = -1.959 Double Precision :: wdnap = 14.98 Double Precision, parameter :: wdpan = -13.83 !Double Precision :: wdndp = 0.0 Double Precision, parameter :: monomer_E = 21000.d0 Double Precision, parameter :: wcm = 1400.0d0 Double Precision, parameter :: Etha = 0.0 !ionicity coefficient! Double Precision :: wGround = -(((monomer_E/wcm)*(1.0-Etha))-((-wdpan-VMU)*Etha)) Double Precision, parameter :: MuAInitial = 1.0 Double Precision :: MuD = 1.0 Double Precision, parameter :: MuA = 1.0 Double Precision :: MuCT = 4.5 !It is along x-axis which is the pi-stack axis! Exciton TDMs are along y exis! Double Precision :: Angle = 0.0 !Angle between dipoles in degree! MuD is considered along y-axis! Double Precision :: Alfadegree = 180.0 !it is a constant to convert from degree to radin! !Double Precision :: Teinter = 2.46 !Double Precision :: Thinter = 2.46 Double Precision :: Tground = 1.6 Double Precision :: Teintra = 0.0 Double Precision :: Thintra = 0.0 Double Precision, parameter :: JDAInitial = 0.4+ Double Precision :: JDA = (MuA/MuAInitial)*JDAInitial Double Precision, parameter :: JDD = 0.0 Double Precision :: JAA = 0.0 Double Precision :: D = 0.0 Double Precision :: wvib = 1.0 Double Precision :: Wmin = 2000.0 !for wavenumber! Double Precision :: Wmax = 49600.0 ! for wavenumber! Double Precision :: Wcut = 19506.590 ! for wavenumber! !Double Precision :: Wmin = 0.0d0 !Double Precision :: Wmax = 15.0d0 Double Precision :: gamLE = 2.26 Double Precision :: gamHE = 2.28 Real, parameter :: PI = 3.1415927 complex*16, parameter :: XJ = (0.d0,1.d0) Double Precision, dimension(:,:), allocatable :: H, HS, OLH1, OLH2, HX, HY Double Precision, dimension(:,:), allocatable :: HSLL, HSLS, HSSL, HSSS Double Precision, dimension(:,:), allocatable :: HLLOff, HLSOff, HSSOff, HSLOff Double Precision, dimension(:,:,:), allocatable :: disorder_elements, wholeP Double Precision, dimension(:,:), allocatable :: H1PE, H2PE, H3PE, HCTnn, HCTnnv Double Precision, dimension(:,:), allocatable :: OL1CT2P, OL22PCT, OL1CT1P, OL21PCT Double Precision, dimension(:,:), allocatable :: OL1CTV2P, OL22PCTV, OL1CTV3P, OL23PCTV Double Precision, dimension(:,:), allocatable :: OL1CTVCT, OL2CTCTV, OL12P1P, OL21P2P Double Precision, dimension(:,:), allocatable :: OL13P2P, OL22P3P Integer, dimension(:,:), allocatable :: indx1 Integer, dimension(:,:,:,:), allocatable :: indx2 Integer, dimension(:,:,:,:), allocatable :: indx3 Integer, dimension(:,:,:,:,:,:), allocatable :: indx4 Integer, dimension(:,:,:,:,:,:), allocatable :: indx5 Double Precision, dimension(:), allocatable :: WORK, OS, Freq, AbX, AbY, Abtot Double Precision, dimension(:,:), allocatable :: LSX, LSY Double Precision, dimension(:), allocatable :: nA, nD1, nD2, SummX, SummY, FX, FY Double Precision, dimension(:), allocatable :: COF1PE, COF2PE, COFCTnn, COFCTnnv, COF3PE Double Precision, dimension(:), allocatable :: COFDpAm, HGround Real(kind=8), dimension(:), allocatable :: Eign !-----------------------------disorderTable-----------------------------------! integer :: Vx, Vy !Double Precision :: rand !Double Precision :: dlarnd !External :: dlarnd integer :: config !integer, parameter :: idist = 3 !integer :: iseed(4)=(/47,3093,1041,77/) !-------------------External Subroutines-----------------------------------------------! EXTERNAL DSYEV EXTERNAL PRINT_MATRIX !EXTERNAL disorder_table !--------------------------------Logical Parameters---------------------------------------! !logical :: LSDiagonal = .true. !logical :: SLDiagonal = .true. !logical :: SSDiagonal = .true. !logical :: LLOffDiagonal = .true. !logical :: LSOffDiagonal = .true. !logical :: SLOffDiagonal = .true. !logical :: SSOffDiagonal = .true. !---------------------Intrinsic Functions-----------------------------------------------! Intrinsic INT, MIN !-----------------------------------------------------------------------------! end module common_variables !---------------------------------------------Start the main program-------------------------------------------! Program Copolymer use common_variables Implicit none !------------------------Start Executive part------------------------------------------! dw = (Wmax - Wmin)/(Z-1) lambP = sqrt(0.5)*lamb !Ground neutral to Excited Cation! lambN = sqrt(0.5)*lamb !!Ground neutral to Excited Anion! lambEP = lamb - lambP !Excited Frenkel to Excited Cation! lambEN = lamb - lambN !Excited Frenkel to Excited Anion! allocate(H(B+1,B+1)) allocate(HX(B,B)) allocate(HY(B,B)) !allocate(disorder_elements(Realization,N,N)) !allocate(wholeP(Realization,N,N)) allocate(OL1CT1P(B3,B1), OL21PCT(B1,B3)) allocate(OL12P1P(B2,B1), OL21P2P(B1,B2)) allocate(OL13P2P(B5,B2), OL22P3P(B2,B5)) allocate(OL1CT2P(B3,B2), OL22PCT(B2,B3)) allocate(OL1CTV2P(B4,B2), OL22PCTV(B2,B4)) allocate(OL1CTV3P(B4,B5), OL23PCTV(B5,B4)) allocate(OL1CTVCT(B4,B3), OL2CTCTV(B3,B4)) !allocate(HSLL(B,B), HSLS(B,B), HSSL(B,B), HSSS(B,B)) !allocate(HLLOff(B,B), HLSOff(B,B), HSSOff(B,B), HSLOff(B,B)) allocate(indx1(N,vibmax+1)) allocate(H1PE(B1,B1), H2PE(B2,B2), HCTnn(B3,B3), HCTnnv(B4,B4), H3PE(B5,B5)) allocate(indx2(N,vibmax,N,vibmax)) allocate(indx3(N,vibmax+1,N,vibmax+1)) allocate(indx4(N,vibmaxT,N,vibmaxT,N,vibmaxT)) allocate(indx5(N,vibmax-1,N,vibmax-1,N,vibmax-1)) allocate(HS(B+1,B+1), Eign(B+1), SummX(B), SummY(B), FX(B), FY(B), Freq(B), OS(B-1), LSX(Z,B), LSY(Z,B)) allocate(nA(Cell), nD1(Cell), nD2(Cell), AbX(Z), AbY(Z), Abtot(Z), HGround(B3)) allocate(COF1PE(B+1), COF2PE(B+1), COFCTnn(B+1), COFCTnnv(B+1), COF3PE(B+1), COFDpAm(B+1)) !-------------------------------HR factor for cation and anion relative to S0 and S1---------------------------! !-----------------------------------------------------------------------------------------------------------------! !-------------------------------Form the disorder table----------------------------------------! !call disorder_table() !---------------------------------Coordinate table for DAD units------------------------------------------------! do i1 = 1,Cell nA(i1) = 3*i1-1 nD1(i1) = 3*i1-2 nD2(i1) = 3*i1 end do print *, 'JDA=',JDA,'' print *, 'wGround=',wGround,'' !---------------------------------------------------------------------------------------------! print *, 'Now Hamiltonian will be formed' !-------------Formation of Hamiltonian----------------------------------------! !----------------------Index-----------------------------------------------------! !--------------------------1PE index------------------------------------! kount1 = 0 do i1 =1,N do j1 =1,(vibmax+1) v1 = j1 - 1 kount1 = kount1 + 1 indx1(i1,j1) = kount1 !print *, 'index <',i1,'/',i2,'> and <',v1,'/',v2,'> is ',indx1(i1,j1,i2,j2),'' end do end do print *, 'kount1=',kount1,'' !-------------------------Formation of the 1PE Hamiltonian--------------------! t1 = 0 do i1 =1,N do j1 =1,(vibmax+1) v1 = j1 - 1 x1 = indx1(i1,j1) do i11 =1,N do j11 =1,(vibmax+1) v11 = j11 - 1 x2 = indx1(i11,j11) if (x2.eq.x1) then if (MOD(i1,2).eq.0) then H1PE(x1,x2) = (wastar + D) + (wvib*(v1*1.0)) else H1PE(x1,x2) = (wdstar + D) + (wvib*(v1*1.0)) end if else if ((iabs(i11-i1)==1).or.(iabs(i11-i1)==N-1).or.(iabs(i11-i1)==2)) then !for Periodic! !else if (iabs(i11-i1)==1) then ! for open! H1PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(t1,v1,lamb)*FC(t1,v11,lamb) else H1PE(x1,x2) = 0.0 end if end do end do end do end do ! print *, 'here is H1PE' !do x1 = 1,B1 ! write(*,59) (H1PE(x1,x2), x2=1,B1) ! end do ! 59 format (8f9.4) !------------------------------2PE index--------------------------------------! kount2 = 0 do i1 =1,N do j1 =1,vibmax v1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle do j2 =1,vibmax v2 = j2 if ((v1+v2)>vibmax) cycle kount2 = kount2 + 1 indx2(i1,j1,i2,j2) = kount2 ! print *, indx2(i1,j1,i2,j2,i3,j3) end do end do end do end do print *, 'kount2=',kount2,'' print *, 'TEST' !------------------------------------Formation of 2PE Hamiltonian----------! !-------------------------------------------------------------------------------! t1 = 0 do i1 =1,N do j1 =1,vibmax v1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle do j2 =1,vibmax v2 = j2 if ((v1+v2)>vibmax) cycle x1 = indx2(i1,j1,i2,j2) do i11 =1,N do j11 =1,vibmax v11 = j11 - 1 do i22=1,N if (i22.eq.i11) cycle do j22 =1,vibmax v22 = j22 if ((v11+v22)>vibmax) cycle x2 = indx2(i11,j11,i22,j22) if (x2.eq.x1) then if (MOD(i1,2).eq.0) then H2PE(x1,x2) = (wastar + D) + (wvib*(v1*1.0 + v2*1.0)) else H2PE(x1,x2) = (wdstar + D) + (wvib*(v1*1.0 + v2*1.0)) end if else if ((iabs(i11-i1)==1).or.(iabs(i11-i1)==N-1).or.(iabs(i11-i1)==2)) then !for periodic! !else if (iabs(i11-i1)==1) then !for open! if ((i22==i2).and.(v22==v2)) then H2PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(t1,v1,lamb)*FC(t1,v11,lamb) else if ((i22==i1).and.(i11==i2)) then H2PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(v22,v1,lamb)*FC(v2,v11,lamb) else H2PE(x1,x2) = 0.0 end if else H2PE(x1,x2) = 0.0 end if end do end do end do end do end do end do end do end do print *, 'TEST2' !print *, 'here is H2PE' !do x1 = 1,B2 !write(*,64) (H2PE(x1,x2), x2=1,B2) !end do !64 format (18f9.4) !---------------------------------------------------------------------------------------------------! !--------------------------------------------3PE index----------------------------------! kount5 = 0 do i1 =1,N do j1 =1,(vibmax-1) v1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle do j2 =1,(vibmax-1) v2 = j2 do i3=1,N !if ((i3.eq.i2).and.(i3.eq.i1)) cycle if ((i3.LE.i2)) cycle if ((i3.eq.i1)) cycle !if (((iabs(i3-i1).ne.1).or.(iabs(i3-i1).ne.N-1)).and.((iabs(i2-i1).ne.1).or.(iabs(i2-i1).ne.N-1))) cycle !if ((iabs(i3-i1).EQ.1).or.(iabs(i3-i1).EQ.N-1).or.(iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then if ((iabs(i3-i1).EQ.1).or.(iabs(i2-i1).EQ.1)) then do j3 =1,(vibmax-1) v3 = j3 if ((v1+v2+v3)>vibmax) cycle kount5 = kount5 + 1 indx5(i1,j1,i2,j2,i3,j3) = kount5 !print *, indx5(i1,j1,i2,j2,i3,j3) !print *, 'index <',i1,'/',i2,'/',i3,'/',v1,'/',v2,'/',v3,'> is ',indx5(i1,j1,i2,j2,i3,j3),'' end do end if end do end do end do end do end do print *, 'kount5=',kount5,'' !---------------------------------------------Formation of 3PE Hamiltonian----------------------------------! t1 = 0 do i1 =1,N do j1 =1,(vibmax-1) v1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle do j2 =1,(vibmax-1) v2 = j2 do i3=1,N if ((i3.LE.i2)) cycle if ((i3.eq.i1)) cycle if ((iabs(i3-i1).EQ.1).or.(iabs(i2-i1).EQ.1)) then do j3 =1,(vibmax-1) v3 = j3 if ((v1+v2+v3)>vibmax) cycle x1 = indx5(i1,j1,i2,j2,i3,j3) do i11 =1,N do j11 =1,(vibmax-1) v11 = j11 - 1 do i22=1,N if (i22.eq.i11) cycle do j22 =1,(vibmax-1) v22 = j22 do i33=1,N if ((i33.LE.i22)) cycle if ((i33.eq.i11)) cycle if ((iabs(i33-i11).EQ.1).or.(iabs(i22-i11).EQ.1)) then do j33 =1,(vibmax-1) v33 = j33 if ((v11+v22+v33)>vibmax) cycle x2 = indx5(i11,j11,i22,j22,i33,j33) if (x2.eq.x1) then if (MOD(i1,2).eq.0) then H3PE(x1,x2) = (wastar + D) + (wvib*(v1*1.0 + v2*1.0 + v3*1.0)) else H3PE(x1,x2) = (wdstar + D) + (wvib*(v1*1.0 + v2*1.0 + v3*1.0)) end if else if ((iabs(i11-i1)==1).or.(iabs(i11-i1)==N-1).or.(iabs(i11-i1)==2)) then !for periodic! !else if (iabs(i11-i1)==1) then !for open! if ((i11==i2).and.(i22==i1).and.(i33==i3).and.(v33==v3)) then H3PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(v2,v11,lamb)*FC(v22,v1,lamb) else if ((i11==i2).and.(i33==i1).and.(i22==i3).and.(v22==v3)) then H3PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(v2,v11,lamb)*FC(v33,v1,lamb) else if ((i11==i3).and.(i22==i1).and.(i33==i2).and.(v33==v2)) then H3PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(v3,v11,lamb)*FC(v22,v1,lamb) else if ((i11==i3).and.(i33==i1).and.(i22==i2).and.(v22==v2)) then H3PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(v3,v11,lamb)*FC(v33,v1,lamb) else if ((i22==i2).and.(v22==v2).and.(i33==i3).and.(v33==v3)) then H3PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(t1,v11,lamb)*FC(t1,v1,lamb) else if ((i22==i3).and.(v22==v3).and.(i33==i2).and.(v33==v2)) then H3PE(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(t1,v11,lamb)*FC(t1,v1,lamb) else H3PE(x1,x2) = 0.0 end if else H3PE(x1,x2) = 0.0 end if end do end if end do end do end do end do end do end do end if end do end do end do end do end do !print *, 'here is H3PE' !do x1 = 1,B5 !write(*,89) (H3PE(x1,x2), x2=1,B5) !end do !89 format (12f9.3) !-----------------------------------------------------------------------------------------------! !------------------------------------------CTnn index-----------------------------------------------! kount3 = 0 do i1=1,N do j1=1,(vibmax+1) vp1 = j1 - 1 do i2=1,N if (i2==i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2=1,(vibmax+1) vn2 = j2 - 1 if ((vp1+vn2)>vibmax) cycle kount3 = kount3 + 1 indx3(i1,j1,i2,j2) = kount3 !print *, indx3(i1,j1,i2,j2) end do end if end do end do end do print *, 'kount3=',kount3,'' !---------------------------------------------Formation of CTnn-----------------------------------! t1 = 0 do i1=1,N do j1=1,(vibmax+1) vp1 = j1 - 1 do i2=1,N if (i2==i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2=1,(vibmax+1) vn2 = j2 - 1 if ((vp1+vn2)>vibmax) cycle x1 = indx3(i1,j1,i2,j2) !print *, 'indexX1 <',i1,'/',vp1,'/',i2,'/',vn2,'>' !print *, 'X1=',x1,'' do i11=1,N do j11=1,(vibmax+1) vp11 = j11 - 1 do i22=1,N if (i22==i11) cycle if ((iabs(i22-i11).EQ.1).or.(iabs(i22-i11).EQ.N-1)) then do j22=1,(vibmax+1) vn22 = j22 - 1 if ((vp11+vn22)>vibmax) cycle x2 = indx3(i11,j11,i22,j22) if (x2.eq.x1) then HCTnn(x1,x2) = (EnergyCT(i1,i2,wdnap,wdpan,N,Sfactor,VMU) + D) + (Wvib*((vp1*1.0)+(vn2*1.0))) else if ((i11==i1).and.(vp11==vp1).and.(iabs((i22-i2))==1)) then !electron moves! HCTnn(x1,x2) = Tefunction(i2,i22,Teintra,N)*FC(t1,vn1,lambN)*FC(t1,vn22,lambN) else if ((i11==i1).and.(vp11==vp1).and.(iabs((i22-i2))==N-1)) then !e-BC! HCTnn(x1,x2) = Tefunction(i2,i22,Teintra,N)*FC(t1,vn1,lambN)*FC(t1,vn22,lambN) else if ((iabs(i11-i1)==1).and.(i22==i2).and.(vn22==vn2)) then !hole moves! HCTnn(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(t1,vp2,lambP)*FC(t1,vp11,lambP) else if ((iabs(i11-i1)==N-1).and.(i22==i2).and.(vn22==vn2)) then !h-BC! HCTnn(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(t1,vp2,lambP)*FC(t1,vp11,lambP) else HCTnn(x1,x2) = 0.0 end if end do end if end do end do end do end do end if end do end do end do !print *, 'here is HCTnn' !do x1 = 1,B3 !write(*,99) (HCTnn(x1,x2), x2=1,B3) !end do !99 format (18f9.3) !----------------------------------------------------------------------------------------------------------------------! !------------------------------Ground state to ionic D+A- electronic coupling---------------------------------------------! !notice that only the CTnn states will have nonezero coupling terms which would be multiplied with FC factors! !CTnnv states will not couple! t1 = 0 do i1=1,N do j1=1,(vibmax+1) vp1 = j1 - 1 do i2=1,N if (i2==i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2=1,(vibmax+1) vn2 = j2 - 1 if ((vp1+vn2)>vibmax) cycle x1 = indx3(i1,j1,i2,j2) if (MOD(i2,2).eq.0) then HGround(x1) = Tground*FC(t1,vp1,lambP)*FC(t1,vn2,lambN) else HGround(x1) = 0.0 end if end do end if end do end do end do !---------------------------------------------------------CTnnv index----------------------------------------------------! kount4 = 0 do i1 =1,N do j1 =1,vibmaxT vp1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2 =1,vibmaxT vn2 = j2 -1 do i3=1,N !if ((i3.eq.i1).or.(i3.eq.i2)) cycle if (i3.EQ.i1) cycle if (i3.EQ.i2) cycle do j3 =1,vibmaxT v3 = j3 if ((vp1+vn2+v3)>vibmaxT) cycle kount4 = kount4 + 1 indx4(i1,j1,i2,j2,i3,j3) = kount4 !print *, 'index <',i1,'/',i2,'/',i3,'/',vp1,'/',vn2,'/',v3,'> is ',indx4(i1,j1,i2,j2,i3,j3),'' end do end do end do end if end do end do end do print *, 'kount4=',kount4,'' !--------------------------------------------Formation of CTnnv----------------------------------! t1 = 0 do i1 =1,N do j1 =1,vibmaxT vp1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2 =1,vibmaxT vn2 = j2 -1 do i3=1,N if (i3.EQ.i1) cycle if (i3.EQ.i2) cycle do j3 =1,vibmaxT v3 = j3 if ((vp1+vn2+v3)>vibmaxT) cycle x1 = indx4(i1,j1,i2,j2,i3,j3) !print *, 'indexX1 <',i1,'/',vp1,'/',i2,'/',vn2,'/',i3,'/',v3,'>' !print *, 'X1=',x1,'' do i11 =1,N do j11 =1,vibmaxT vp11 = j11 - 1 do i22=1,N if (i22.eq.i11) cycle if ((iabs(i22-i11).EQ.1).or.(iabs(i22-i11).EQ.N-1)) then do j22 =1,vibmaxT vn22 = j22 - 1 do i33=1,N if (i33.EQ.i11) cycle if (i33.EQ.i22) cycle do j33 =1,vibmaxT v33 = j33 if ((vp11+vn22+v33)>vibmaxT) cycle x2 = indx4(i11,j11,i22,j22,i33,j33) !print *, 'indexX2 <',i11,'/',vp11,'/',i22,'/',vn22,'/',i33,'/',v33,'>' !print *, 'X1=',x2,'' if (x2.eq.x1) then HCTnnv(x1,x2) = (EnergyCT(i1,i2,wdnap,wdpan,N,Sfactor,VMU) + D) + (Wvib*((vp1*1.0)+(vn2*1.0)+(v3*1.0))) else if ((i11==i1).and.(vp11==vp1).and.(iabs((i22-i2))==1)) then !electron moves! if ((i33==i3).and.(v33==v3)) then HCTnnv(x1,x2) = Tefunction(i2,i22,Teintra,N)*FC(t1,vn1,lambN)*FC(t1,vn22,lambN) else if ((i33==i2).and.(i22==i3)) then HCTnnv(x1,x2) = Tefunction(i2,i22,Teintra,N)*FC(v33,vn1,lambN)*FC(v3,vn22,lambN) end if else if ((i11==i1).and.(vp11==vp1).and.(iabs((i22-i2))==N-1)) then !electron BC! if ((i33==i3).and.(v33==v3)) then HCTnnv(x1,x2) = Tefunction(i2,i22,Teintra,N)*FC(t1,vn1,lambN)*FC(t1,vn22,lambN) else if ((i33==i2).and.(i22==i3)) then HCTnnv(x1,x2) = Tefunction(i2,i22,Teintra,N)*FC(v33,vn1,lambN)*FC(v3,vn22,lambN) end if else if ((iabs(i11-i1)==1).and.(i22==i2).and.(vn22==vn2)) then !hole moves! if ((i33==i3).and.(v33==v3)) then HCTnnv(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(t1,vp1,lambP)*FC(t1,vp11,lambP) else if ((i33==i1).and.(i11==i3)) then HCTnnv(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(v33,vp1,lambP)*FC(v3,vp11,lambP) end if else if ((iabs(i11-i1)==N-1).and.(i22==i2).and.(vn22==vn2)) then !hole BC! if ((i33==i3).and.(v33==v3)) then HCTnnv(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(t1,vp1,lambP)*FC(t1,vp11,lambP) else if ((i33==i1).and.(i11==i3)) then HCTnnv(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(v33,vp1,lambP)*FC(v3,vp11,lambP) end if else HCTnnv(x1,x2) = 0.0 end if end do end do end do end if end do end do end do end do end do end do end if end do end do end do ! print *, 'here is HCTnnv' ! do x1 = 1,B4 ! write(*,65) (HCTnnv(x1,x2), x2=1,B4) ! end do ! 65 format (4f9.3) !-----------------------------------------------------------------------------------------------! !----------------------------------------------OL12P1P and OL21P2P--------------------------------! t1 = 0 do i1 =1,N do j1 =1,vibmax v1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle do j2 =1,vibmax v2 = j2 if ((v1+v2)>vibmax) cycle x1 = indx2(i1,j1,i2,j2) do i11 =1,N do j11 =1,(vibmax+1) v11 = j11 - 1 x2 = indx1(i11,j11) if ((iabs(i11-i1)==1).or.(iabs(i11-i1)==N-1).or.(iabs(i11-i1)==2)) then if (i11==i2) then OL12P1P(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(t1,v1,lamb)*FC(v2,v11,lamb) OL21P2P(x2,x1) = OL12P1P(x1,x2) else OL12P1P(x1,x2) = 0.0 OL21P2P(x2,x1) = OL12P1P(x1,x2) end if end if end do end do end do end do end do end do !print *, 'here is OL12P1P' !do x1 = 1,B2 !write(*,76) (OL12P1P(x1,x2), x2=1,B1) !end do !76 format (9f9.3) !print *, 'here is OL21P2P' !do x1 = 1,B1 !write(*,46) (OL21P2P(x1,x2), x2=1,B2) !end do !46 format (18f9.3) !----------------------------------------------OL13P2P and OL22P3P--------------------------------! t1 = 0 do i1 =1,N do j1 =1,(vibmax-1) v1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle do j2 =1,(vibmax-1) v2 = j2 do i3=1,N if ((i3.LE.i2)) cycle if ((i3.eq.i1)) cycle if ((iabs(i3-i1).EQ.1).or.(iabs(i2-i1).EQ.1)) then do j3 =1,(vibmax-1) v3 = j3 if ((v1+v2+v3)>vibmax) cycle x1 = indx5(i1,j1,i2,j2,i3,j3) do i11 =1,N do j11 =1,vibmax v11 = j11 - 1 do i22=1,N if (i22.eq.i11) cycle do j22 =1,vibmax v22 = j22 if ((v11+v22)>vibmax) cycle x2 = indx2(i11,j11,i22,j22) if ((iabs(i11-i1)==1).or.(iabs(i11-i1)==N-1).or.(iabs(i11-i1)==2)) then if ((i11==i2).and.(i22==i3).and.(v22==v3)) then OL13P2P(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(t1,v1,lamb)*FC(v2,v11,lamb) OL22P3P(x2,x1) = OL13P2P(x1,x2) else if ((i11==i3).and.(i22==i2).and.(v22==v2)) then OL13P2P(x1,x2) = Jfunction(i1,i11,JDA,JDD,JAA,N)*FC(t1,v1,lamb)*FC(v3,v11,lamb) OL22P3P(x2,x1) = OL13P2P(x1,x2) else OL13P2P(x1,x2) = 0.0 OL22P3P(x2,x1) = OL13P2P(x1,x2) end if end if end do end do end do end do end do end if end do end do end do end do end do !print *, 'here is OL13P2P' !do x1 = 1,B5 !write(*,43) (OL13P2P(x1,x2), x2=1,B2) !end do !43 format (9f9.3) !print *, 'here is OL21P2P' !do x1 = 1,B2 !write(*,38) (OL22P3P(x1,x2), x2=1,B5) !end do !38 format (18f9.3) !----------------------------------------------OL1CT1P and OL21PCT--------------------------------! t1 = 0 do i1=1,N do j1=1,(vibmax+1) vp1 = j1 - 1 do i2=1,N if (i2==i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2=1,(vibmax+1) vn2 = j2 - 1 if ((vp1+vn2)>vibmax) cycle x1 = indx3(i1,j1,i2,j2) do i11=1,N do j11=1,(vibmax+1) x2 = indx1(i11,j11) v11 = j11 - 1 if ((i11==i1).and.(iabs((i11-i2))==1)) then !e-moves! OL1CT1P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(t1,vn2,lambN) OL21PCT(x2,x1) = OL1CT1P(x1,x2) else if ((i11==i1).and.(iabs((i11-i2))==N-1)) then !e-BC! OL1CT1P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(t1,vn2,lambN) OL21PCT(x2,x1) = OL1CT1P(x1,x2) else if ((iabs(i11-i1)==1).and.(i11==i2)) then !h-moves! OL1CT1P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(t1,vp1,lambP) OL21PCT(x2,x1) = OL1CT1P(x1,x2) else if ((iabs(i11-i1)==N-1).and.(i11==i2)) then !h-BC! OL1CT1P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(t1,vp1,lambP) OL21PCT(x2,x1) = OL1CT1P(x1,x2) else OL1CT1P(x1,x2) = 0.0 OL21PCT(x2,x1) = OL1CT1P(x1,x2) end if end do end do end do end if end do end do end do !print *, 'here is OL1CT1P' !do x1 = 1,B3 !write(*,94) (OL1CT1P(x1,x2), x2=1,B1) !end do !94 format (9f9.3) !print *, 'here is OL21PCT' !do x1 = 1,B1 ! write(*,27) (OL21PCT(x1,x2), x2=1,B3) !end do ! 27 format (16f9.3) !--------------------------------------------------------OL1CT2P and OL22PCT--------------------------! t1 = 0 do i1=1,N do j1=1,(vibmax+1) vp1 = j1 - 1 do i2=1,N if (i2==i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2=1,(vibmax+1) vn2 = j2 - 1 if ((vp1+vn2)>vibmax) cycle x1 = indx3(i1,j1,i2,j2) do i11=1,N do j11=1,vibmax v11 = j11 - 1 do i22=1,N if (i22==i11) cycle do j22=1,vibmax v22 = j22 if ((v11+v22)>vibmax) cycle x2 = indx2(i11,j11,i22,j22) if ((i11==i1).and.(iabs(i11-i2)==1).and.(i22==i2)) then !e-moves! OL1CT2P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(v22,vn2,lambN) OL22PCT(x2,x1) = OL1CT2P(x1,x2) else if ((i11==i1).and.(iabs(i11-i2)==N-1).and.(i22==i2)) then !e-BC! OL1CT2P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(v22,vn2,lambN) OL22PCT(x2,x1) = OL1CT2P(x1,x2) else if ((iabs(i11-i1)==1).and.(i11==i2).and.(i22==i1)) then !h-moves! OL1CT2P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(v22,vp1,lambP) OL22PCT(x2,x1) = OL1CT2P(x1,x2) else if ((iabs(i11-i1)==N-1).and.(i11==i2).and.(i22==i1)) then !h-BC! OL1CT2P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(v22,vp1,lambP) OL22PCT(x2,x1) = OL1CT2P(x1,x2) else OL1CT2P(x1,x2) = 0.0 OL22PCT(x2,x1) = OL1CT2P(x1,x2) end if end do end do end do end do end do end if end do end do end do !print *, 'here is OL1CT2P' !do x1 = 1,B3 !write(*,83) (OL1CT2P(x1,x2), x2=1,B2) !end do !83 format (6f9.3) !print *, 'here is OL22PCT' !do x1 = 1,B2 ! write(*,46) (OL22PCT(x1,x2), x2=1,B3) !end do !46 format (12f9.3) !-------------------------------------------------OL1CTV2P and OL22PCTV-----------------------------! t1 = 0 do i1 =1,N do j1 =1,vibmaxT vp1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2 =1,vibmaxT vn2 = j2 -1 do i3=1,N if (i3.EQ.i1) cycle if (i3.EQ.i2) cycle do j3 =1,vibmaxT v3 = j3 if ((vp1+vn2+v3)>vibmaxT) cycle x1 = indx4(i1,j1,i2,j2,i3,j3) do i11=1,N do j11=1,vibmax v11 = j11 - 1 do i22=1,N if (i22==i11) cycle do j22=1,vibmax v22 = j22 if ((v11+v22)>vibmax) cycle x2 = indx2(i11,j11,i22,j22) if ((i11==i1).and.(iabs(i1-i2)==1).and.(i22==i3).and.(v22==v3)) then !e-moves! OL1CTV2P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(t1,vn2,lambN) OL22PCTV(x2,x1) = OL1CTV2P(x1,x2) else if ((i11==i1).and.(iabs(i1-i2)==N-1).and.(i22==i3).and.(v22==v3)) then !e-BC! OL1CTV2P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(t1,vn2,lambN) OL22PCTV(x2,x1) = OL1CTV2P(x1,x2) else if ((iabs(i1-i2)==1).and.(i11==i2).and.(i22==i3).and.(v22==v3)) then !h-moves! OL1CTV2P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(t1,vp1,lambP) OL22PCTV(x2,x1) = OL1CTV2P(x1,x2) else if ((iabs(i1-i2)==N-1).and.(i11==i2).and.(i22==i3).and.(v22==v3)) then !h-BC! OL1CTV2P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(t1,vp1,lambP) OL22PCTV(x2,x1) = OL1CTV2P(x1,x2) else OL1CTV2P(x1,x2) = 0.0 OL22PCTV(x2,x1) = OL1CTV2P(x1,x2) end if end do end do end do end do end do end do end do end if end do end do end do !print *, 'here is OL1CTV2P' !do x1 = 1,B4 !write(*,43) (OL1CTV2P(x1,x2), x2=1,B2) !end do !43 format (6f9.3) !print *, 'here is OL22PCTV' !do x1 = 1,B2 !write(*,95) (OL22PCTV(x1,x2), x2=1,B4) !end do !95 format (4f9.3) !----------------------------------------------------OL1CTV3P and OL23PCTV------------------------------------------! t1 = 0 do i1 =1,N do j1 =1,vibmaxT vp1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2 =1,vibmaxT vn2 = j2 -1 do i3=1,N if (i3.EQ.i1) cycle if (i3.EQ.i2) cycle do j3 =1,vibmaxT v3 = j3 if ((vp1+vn2+v3)>vibmaxT) cycle x1 = indx4(i1,j1,i2,j2,i3,j3) do i11 =1,N do j11 =1,(vibmax-1) v11 = j11 - 1 do i22=1,N if (i22.eq.i11) cycle do j22 =1,(vibmax-1) v22 = j22 do i33=1,N if ((i33.LE.i22)) cycle if ((i33.eq.i11)) cycle if ((iabs(i33-i11).EQ.1).or.(iabs(i22-i11).EQ.1)) then do j33 =1,(vibmax-1) v33 = j33 if ((v11+v22+v33)>vibmax) cycle x2 = indx5(i11,j11,i22,j22,i33,j33) if ((i11==i1).and.(iabs((i2-i11))==1)) then !electron moves! if ((i22==i2).and.(i33==i3).and.(v33==v3)) then OL1CTV3P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(v22,vn2,lambN) OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) else if ((i22==i3).and.(i33==i2).and.(v22==v3)) then OL1CTV3P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(v33,vn2,lambN) OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) end if else if ((i11==i1).and.(iabs((i2-i11))==N-1)) then !electron BC! if ((i22==i2).and.(i33==i3).and.(v33==v3)) then OL1CTV3P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(v22,vn2,lambN) OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) else if ((i22==i3).and.(i33==i2).and.(v22==v3)) then OL1CTV3P(x1,x2) = Tefunction(i2,i11,Teintra,N)*FC(vp1,v11,lambEP)*FC(v33,vn2,lambN) OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) end if else if ((iabs(i11-i1)==1).and.(i11==i2)) then !hole moves! if ((i22==i1).and.(i33==i3).and.(v33==v3)) then OL1CTV3P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(v22,vp1,lambP) OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) else if ((i33==i1).and.(i22==i3).and.(v22==v3)) then OL1CTV3P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(v33,vp1,lambP) OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) end if else if ((iabs(i11-i1)==N-1).and.(i11==i2)) then !hole BC! if ((i22==i1).and.(i33==i3).and.(v33==v3)) then OL1CTV3P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(v22,vp1,lambP) OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) else if ((i33==i1).and.(i22==i3).and.(v22==v3)) then OL1CTV3P(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(vn2,v11,lambEN)*FC(v33,vp1,lambP) OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) end if else OL1CTV3P(x1,x2) = 0.0 OL23PCTV(x2,x1) = OL1CTV3P(x1,x2) end if end do end if end do end do end do end do end do end do end do end do end if end do end do end do !print *, 'here is OL1CTV3P' !do x1 = 1,B4 !write(*,92) (OL1CTV2P(x1,x2), x2=1,B5) !end do !92 format (3f9.3) !print *, 'here is OL23PCTV' !do x1 = 1,B5 !write(*,83) (OL22PCTV(x1,x2), x2=1,B4) !end do !83 format (16f9.3) !----------------------------------------------OL1CTVCT and OL2CTCTV---------------------------------------------------! t1 = 0 do i1 =1,N do j1 =1,vibmaxT vp1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2 =1,vibmaxT vn2 = j2 -1 do i3=1,N if (i3.EQ.i1) cycle if (i3.EQ.i2) cycle do j3 =1,vibmaxT v3 = j3 if ((vp1+vn2+v3)>vibmaxT) cycle x1 = indx4(i1,j1,i2,j2,i3,j3) do i11=1,N do j11=1,(vibmax+1) vp11 = j11 - 1 do i22=1,N if (i22==i11) cycle if ((iabs(i22-i11).EQ.1).or.(iabs(i22-i11).EQ.N-1)) then do j22=1,(vibmax+1) vn22 = j22 - 1 if ((vp11+vn22)>vibmax) cycle x2 = indx3(i11,j11,i22,j22) if ((i11==i1).and.(vp11==vp1).and.(iabs((i2-i22))==1).and.(i22==i3)) then !electron moves! OL1CTVCT(x1,x2) = Tefunction(i2,i22,Teintra,N)*FC(v3,vn22,lambN)*FC(t1,vn1,lambN) OL2CTCTV(x2,x1) = OL1CTVCT(x1,x2) else if ((i11==i1).and.(vp11==vp1).and.(iabs((i2-i22))==N-1).and.(i22==i3)) then !electron BC! OL1CTVCT(x1,x2) = Tefunction(i2,i22,Teintra,N)*FC(v3,vn22,lambN)*FC(t1,vn1,lambN) OL2CTCTV(x2,x1) = OL1CTVCT(x1,x2) else if ((iabs(i11-i1)==1).and.(i22==i2).and.(vn22==vn1).and.(i11==i3)) then !hole moves! OL1CTVCT(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(v3,vp11,lambP)*FC(t1,vp1,lambP) OL2CTCTV(x2,x1) = OL1CTVCT(x1,x2) else if ((iabs(i11-i1)==N-1).and.(i22==i2).and.(vn22==vn1).and.(i11==i3)) then !hole BC! OL1CTVCT(x1,x2) = Thfunction(i1,i11,Thintra,N)*FC(v3,vp11,lambP)*FC(t1,vp1,lambP) OL2CTCTV(x2,x1) = OL1CTVCT(x1,x2) else OL1CTVCT(x1,x2) = 0.0 OL2CTCTV(x2,x1) = OL1CTVCT(x1,x2) end if end do end if end do end do end do end do end do end do end if end do end do end do print *, 'TEST3' !print *, 'here is OL1CTVCT' !do x1 = 1,B4 !write(*,18) (OL1CTVCT(x1,x2), x2=1,B3) !end do !18 format (12f9.3) !print *, 'here is OL2CTCTV' !do x1 = 1,B3 !write(*,16) (OL2CTCTV(x1,x2), x2=1,B4) !end do !16 format (4f9.3) !----------------------------Now, Form the total Hamiltonian---------------------------------! do i =1,(B+1) do j =1,(B+1) if (i==1) then if ((j>(B1+B2+1)).and.(j<=(B1+B2+B3+1))) then x2 = j - (B1+B2+1) H(i,j) = Hground(x2) H(j,i) = H(i,j) else if (j==1) then H(i,j) = wGround else H(i,j) = 0.0 H(j,i) = 0.0 end if end if end do end do do i =1,B do j =1,B if (i<=B1) then if (j<=B1) then H(i+1,j+1) = H1PE(i,j) else if ((j>B1).and.(j<=(B1+B2))) then x2 = j - B1 H(i+1,j+1) = OL21P2P(i,x2) else if ((j>(B1+B2)).and.(j<=(B1+B2+B3))) then x2 = j - (B1+B2) H(i+1,j+1) = OL21PCT(i,x2) else if ((j>(B1+B2+B3)).and.(j<=(B1+B2+B3+B4))) then x2 = j - (B1+B2+B3) H(i+1,j+1) = 0.0 else if ((j>(B1+B2+B3+B4)).and.(j<=B)) then x2 = j - (B1+B2+B3+B4) H(i+1,j+1) = 0.0 !print *, 'TEST3' end if !print *, 'TEST4' else if ((i>B1).and.(i<=(B1+B2))) then x1 = i - B1 if (j<=B1) then H(i+1,j+1) = OL12P1P(x1,j) !print *, 'TEST5' else if ((j>B1).and.(j<=(B1+B2))) then x2 = j - B1 H(i+1,j+1) = H2PE(x1,x2) !print *, 'TEST6' else if ((j>(B1+B2)).and.(j<=(B1+B2+B3))) then x2 = j - (B1+B2) H(i+1,j+1) = OL22PCT(x1,x2) !print *, 'TEST7' else if ((j>(B1+B2+B3)).and.(j<=(B1+B2+B3+B4))) then x2 = j - (B1+B2+B3) H(i+1,j+1) = OL22PCTV(x1,x2) !print *, 'TEST8' else if ((j>(B1+B2+B3+B4)).and.(j<=B)) then x2 = j - (B1+B2+B3+B4) H(i+1,j+1) = OL22P3P(x1,x2) !H(i,j) = 0.0 !for test! !print *, 'TEST9' end if else if ((i>(B1+B2)).and.(i<=(B1+B2+B3))) then x1 = i - (B1+B2) if (j<=B1) then H(i+1,j+1) = OL1CT1P(x1,j) !print *, 'TEST10' else if ((j>B1).and.(j<=(B1+B2))) then x2 = j - B1 H(i+1,j+1) = OL1CT2P(x1,x2) !print *, 'TEST11' else if ((j>(B1+B2)).and.(j<=(B1+B2+B3))) then x2 = j - (B1+B2) H(i+1,j+1) = HCTnn(x1,x2) !print *, 'TEST12' else if ((j>(B1+B2+B3)).and.(j<=(B1+B2+B3+B4))) then x2 = j - (B1+B2+B3) H(i+1,j+1) = OL2CTCTV(x1,x2) !print *, 'TEST8' else if ((j>(B1+B2+B3+B4)).and.(j<=B)) then x2 = j - (B1+B2+B3+B4) H(i+1,j+1) = 0.0 !print *, 'TEST13' end if !print *, 'TEST6' else if ((i>(B1+B2+B3)).and.(i<=(B1+B2+B3+B4))) then x1 = i - (B1+B2+B3) if (j<=B1) then H(i+1,j+1) = 0.0 else if ((j>B1).and.(j<=(B1+B2))) then x2 = j - B1 H(i+1,j+1) = OL1CTV2P(x1,x2) else if ((j>(B1+B2)).and.(j<=(B1+B2+B3))) then x2 = j - (B1+B2) H(i+1,j+1) = OL1CTVCT(x1,x2) else if ((j>(B1+B2+B3)).and.(j<=(B1+B2+B3+B4))) then x2 = j - (B1+B2+B3) H(i+1,j+1) = HCTnnv(x1,x2) else if ((j>(B1+B2+B3+B4)).and.(j<=B)) then x2 = j - (B1+B2+B3+B4) H(i+1,j+1) = OL1CTV3P(x1,x2) end if else if ((i>(B1+B2+B3+B4)).and.(i<=B)) then x1 = i - (B1+B2+B3+B4) if (j<=B1) then H(i+1,j+1) = 0.0 else if ((j>B1).and.(j<=(B1+B2))) then x2 = j - B1 H(i+1,j+1) = OL13P2P(x1,x2) !H(i,j) = 0.0 !for test! else if ((j>(B1+B2)).and.(j<=(B1+B2+B3))) then x2 = j - (B1+B2) H(i+1,j+1) = 0.0 else if ((j>(B1+B2+B3)).and.(j<=(B1+B2+B3+B4))) then x2 = j - (B1+B2+B3) H(i+1,j+1) = OL23PCTV(x1,x2) else if ((j>(B1+B2+B3+B4)).and.(j<=B)) then x2 = j - (B1+B2+B3+B4) H(i+1,j+1) = H3PE(x1,x2) end if end if end do end do !print *, 'here is the total H' !do x1 = 1,(B+1) !write(*,39) (H(x1,x2), x2=1,(B+1)) !end do !39 format (13f5.1) !-------------------------------Check each part of Hamiltonian-------------------------------! !print *, 'kount1=',kount1,'' !print *, 'kount2=',kount2,'' !print *, 'kount3=',kount3,'' !print *, 'kount4=',kount4,'' !print *, 'kount5=',kount5,'' print *, 'Sum of kounts=',kount1+kount2+kount3+kount4+kount5,'' print *, 'Dimension of Hamiltonian=',B+1,'' !-----------------------------Check the Hamiltonian-----------------------------------------! !print *, 'here is the total H' ! do x1 = 1,B ! write(*,1) (H(x1,x2), x2=1,B) ! end do ! 1 format (12f5.1) !-------------------------------Is it correct Hamiltonian?------------------------------------! !print *, 'Is it correct Hamiltonian? if yes, type Y, and if no, type N' !read *, Answer ! if (Answer /= 'Y') then ! continue ! else ! stop ! end if !------------------------------Copy the original Hamiltonian--------------------------------------! do x1=1,(B+1) do x2=1,(B+1) HS(x1,x2) = H(x1,x2) !HSLL(x1,x2) = H(x1,x2) !HSLS(x1,x2) = H(x1,x2) !HSSL(x1,x2) = H(x1,x2) !HSSS(x1,x2) = H(x1,x2) !HLLOff(x1,x2) = H(x1,x2) !HSLOff(x1,x2) = H(x1,x2) !HLSOff(x1,x2) = H(x1,x2) !HSSOff(x1,x2) = H(x1,x2) end do end do !--------------------------------Let's add Disorder-------------------------------------------! !Print *, "Diagonal Disorder:)" !call LLDiagonal() !call SLDiagonal() !call LSDiagonal() !call SSDiagonal() !print *, "just finished diagonal disorder:))" !print *, "Off-diagonal disorder--Ugly Face" !call LLOffDiagonal() !call SLOffDiagonal() !call LSOffDiagonal() !call SSOffDiagonal() !print *, "just finished Off-diagonal disorder:))" !--------------------------------------------------------------------------------------------! !-----------------------------NOW, Diagonalize the Hamiltonian------------------------------! Write(*,*) 'DSYEV program results' !-----------------------------------------------------------! LWORK = -1 allocate(WORK(3*(B)-1)) Call DSYEV('Vector', 'Upper', (B+1), H, (B+1), Eign, WORK, LWORK, INFO) LWORK=WORK(1) deallocate(WORK) allocate(WORK(LWORK)) Call DSYEV('Vector', 'Upper', (B+1), H, (B+1), Eign, WORK, LWORK, INFO) !--------------------------------------------------------------------------! !----------------------------Check for Convergence--------------------------------------! if (INFO.gt.0) then write(*,*) 'The algorithm failed to compute the eigenvalue' stop end if !---------------------------Print the eigenvalues and eigenvectors---------------------------! !CALL PRINT_MATRIX('eigenvalues', 1, B, Eign, 1) !CALL PRINT_MATRIX('Eigenvectores--stored columnwise', B, B, H, LDA) !------------------------------------So, the results are-----------------------------------------! !-----------------------index 1 and 2 goes over x and y direction 1:y and 2:x----------------------! do x4 = 1,13 print *, Eign(x4) end do !-----------------------------------------------CT Character-----------------------------------------! do x2 = 1,B COF1PE (x2) = 0.0 COF2PE(x2) = 0.0 COFCTnn(x2) = 0.0 COFCTnnv(x2) = 0.0 COF3PE(x2) = 0.0 do x4 = 1, B1 COF1PE(x2) = COF1PE(x2) + H(x4,x2)*H(x4,x2) end do do x4 = 1, B2 v4 = x4 + B1 COF2PE(x2) = COF2PE(x2) + H(v4,x2)*H(v4,x2) end do do x4 = 1, B3 v4 = x4 + (B1+B2) COFCTnn(x2) = COFCTnn(x2) + H(v4,x2)*H(v4,x2) end do do x4 = 1, B4 v4 = x4 + (B1+B2+B3) COFCTnnv(x2) = COFCTnnv(x2) + H(v4,x2)*H(v4,x2) end do do x4 = 1, B5 v4 = x4 + (B1+B2+B3+B4) COF3PE(x2) = COF3PE(x2) + H(v4,x2)*H(v4,x2) end do end do !----------------------------------------the probability of D+A-; among all the CT states---------------------------------------! do X2 = 1, B COFDpAm (x2) = 0.0 !Go over CTnn! do i1=1,N do j1=1,(vibmax+1) vp1 = j1 - 1 do i2=1,N if (i2==i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2=1,(vibmax+1) vn2 = j2 - 1 if ((vp1+vn2)>vibmax) cycle x1 = indx3(i1,j1,i2,j2) x4 = (B1+B2) + x1 if (MOD(i2,2).eq.0) then COFDpAm(x2) = COFDpAm(x2) + H(x4,x2)*H(x4,x2) end if end do end if end do end do end do ! Go over CTnnv! do i1 =1,N do j1 =1,vibmaxT vp1 = j1 - 1 do i2=1,N if (i2.eq.i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2 =1,vibmaxT vn2 = j2 -1 do i3=1,N if (i3.EQ.i1) cycle if (i3.EQ.i2) cycle do j3 =1,vibmaxT v3 = j3 if ((vp1+vn2+v3)>vibmaxT) cycle x1 = indx4(i1,j1,i2,j2,i3,j3) x4 = (B1+B2+B3) + x1 if (MOD(i2,2).eq.0) then COFDpAm(x2) = COFDpAm(x2) + H(x4,x2)*H(x4,x2) end if end do end do end do end if end do end do end do end do !-------------------------------------------------------------------------------------------------------! print *, 'now for OSillator strength' do x2 =2,(B+1) x3 = x2 - 1 !Go over 1PE states! do i1 = 1,N do j1 =1,(vibmax+1) x1 = indx1(i1,j1) x5 = x1 x1 = x1 + 1 v1 = j1 - 1 t1 = 0 if (MOD(i1,2).eq.0) then x4 = 2 ! HX(x5,x3) = Dipoles(x4,i1,MuA,MuD,Angle,Alfadegree)*((H(1,1)*H(x1,x2)*FC(t1,v1,lamb))+(H(x1,1)*H(1,x2)*FC(t1,v1,lamb)) & ! +(H(x1,1)*H(x1,x2))) HX(x5,x3) = 0.0 x4 = 1 HY(x5,x3) = Dipoles(x4,i1,MuA,MuD,Angle,Alfadegree)*((H(1,1)*H(x1,x2)*FC(t1,v1,lamb))+(H(x1,1)*H(1,x2)*FC(t1,v1,lamb)) & +(H(x1,1)*H(x1,x2))) !print *, 'FCA and v1 <',FC(t1,v1,lamb),'/',v1,'>' else x4 = 2 ! HX(x5,x3) = Dipoles(x4,i1,MuA,MuD,Angle,Alfadegree)*((H(1,1)*H(x1,x2)*FC(t1,v1,lamb))+(H(x1,1)*H(1,x2)*FC(t1,v1,lamb)) & ! +(H(x1,1)*H(x1,x2))) HX(x5,x3) = 0.0 x4 = 1 HY(x5,x3) = Dipoles(x4,i1,MuA,MuD,Angle,Alfadegree)*((H(1,1)*H(x1,x2)*FC(t1,v1,lamb))+(H(x1,1)*H(1,x2)*FC(t1,v1,lamb)) & +(H(x1,1)*H(x1,x2))) !print *, 'FCD and v1 <',FC(t1,v1,lamb),'/',v1,'>' end if end do end do !Go over CTnn states! do i1=1,N do j1=1,(vibmax+1) vp1 = j1 - 1 do i2=1,N if (i2==i1) cycle if ((iabs(i2-i1).EQ.1).or.(iabs(i2-i1).EQ.N-1)) then do j2=1,(vibmax+1) vn2 = j2 - 1 if ((vp1+vn2)>vibmax) cycle x1 = indx3(i1,j1,i2,j2) x5 = x1 + B1 !x1 = x1 + B1 x1 = (x1 + 1 + B1 + B2) if (MOD(i2,2).eq.0) then HX(x5,x3) = MuCT*H(x1,1)*H(x1,x2) !HY(x5,x3) = MuCT*H(x1,1)*H(x1,x2) else HX(x5,x3) = 0.0 HY(x5,x3) = 0.0 end if end do end if end do end do end do end do !................To calculate Oscillator Strength........................................................! do x2=1,B SummX (x2) = 0.0 SummY (x2) = 0.0 do x1=1,B SummX(x2) = SummX(x2) + HX(x1,x2) SummY(x2) = SummY(x2) + HY(x1,x2) end do end do !--------------------------------------! do x2=1,B SummX(x2) = SummX(x2)*SummX(x2) SummY(x2) = SummY(x2)*SummY(x2) end do !---------------------------------------! print *, 'the values of oscillator strength are as follows: ' do x2=1,B FX(x2) = SummX(x2) FY(x2) = SummY(x2) !print *, 'os and f <',SummX(x2),'/',SummY(x2),'>' !Freq(x2) = Eign(x2)*wcm*1.d0 + monomer_E*1.0 !Freq(x2) = Eign(x2)*wcm*1.d0 !F(x2) = Freq(x2)*F(x2) x3 = x2 + 1 Freq(x2) = (Eign(x3)-Eign(1))*wcm*1.d0 !for wavenumber! !print *, 'os and f <',SummX(x2),'/',SummY(x2),'/',Freq(x2),'>' !Freq(x2) = Eign(x2)*wcm*1.d0/8065.0 + monomer_E*1.0 FX(x2) = Freq(x2)*FX(x2) !here I multiply with frequency! FY(x2) = Freq(x2)*FY(x2) !print *, Freq(x2) ! print *, 'freq and f <',Freq(x2),'/',F(x2),'>' !F(x2) = Freq(x2)*F(x2) !write(*,8) F(x2) end do ! 8 format (1f5.1) !------------------------------------------------------------------------------------------! !-----------------------------------Form the Line Shape Matrix, Absorption---------------------! !Fmax = MAXVAL(FY) !do x4 = 1,B !print *, 'Wavelength and w and f <',(10000000.0/Freq(x4)),'/',(x4*1.0),'/',(FY(x4)/Fmax),'>' !end do gamLE = gamLE*wcm gamHE = gamHE*wcm do x1=1,Z w = Wmin + ((x1-1)*dw) !w = w*wcm*1.d0 + monomer_E*1.0 !w = w*wcm*1.d0 AbX(x1) = 0.0 AbY(x1) = 0.0 do x2=1,B !LS(x1,x2) = DEXP(-((w-Eign(x2))**2)/(gam**2)) if (Freq(x2).LE.Wcut) then LSX(x1,x2) = (1.0d0/gamLE)*DEXP(-((w-Freq(x2))**2)/(gamLE**2)) LSY(x1,x2) = (1.0d0/gamLE)*DEXP(-((w-Freq(x2))**2)/(gamLE**2)) else LSX(x1,x2) = (1.0d0/gamHE)*DEXP(-((w-Freq(x2))**2)/(gamHE**2)) LSY(x1,x2) = (1.0d0/gamHE)*DEXP(-((w-Freq(x2))**2)/(gamHE**2)) end if !print *, 'freq and w and f <',Freq(x2),'/',w,'/',F(x2),'>' !LS(x1,x2) = (gam**2)/(((w-Eign(x2))**2)+(gam**2)) !F(x2) = F(x2)*abs((w-Eign(x2))) !Ab(x1) = Ab(x1) + (F(x2)*(Eign(x2))*LS(x1,x2)) !Ab(x1) = Ab(x1) + (F(x2)*(Freq(x2) - monomer_E*1.0)*LS(x1,x2)) AbX(x1) = AbX(x1) + (FX(x2)*LSX(x1,x2)) AbY(x1) = AbY(x1) + (FY(x2)*LSY(x1,x2)) !print *, 'freq and w and f <',Freq(x2),'/',w,'/',F(x2),'/',Ab(x1),'>' end do AbX(x1) = (1.0/N)*AbX(x1) AbY(x1) = (1.0/N)*AbY(x1) Abtot(x1) = AbX(x1) + AbY(x1) !print *, 'w and abs <',w,'/',Ab(x1),'>' ! print *, Ab(x1) end do !renormalize spectrum! AbmaxX = MAXVAL(AbX) AbmaxY = MAXVAL(AbY) Abmaxtot = MAXVAL(Abtot) do x1=1,Z AbX(x1) = AbX(x1)/AbmaxX AbY(x1) = AbY(x1)/AbmaxY Abtot(x1) = Abtot(x1)/Abmaxtot end do !print *, 'check!' !do x1= 1,B !print *, Eign(x1) !end do !-----------------------------------------Show the absorption for each frequency-------------------! print *, 'Here is the spectrum' open(UNIT=4,FILE="COPVKSVMU1.txt",STATUS="OLD",ACTION="READWRITE") do x1=1,Z w = Wmin + ((x1-1)*(dw)) weV = w/8065.0 !weV = w*1400.0/8065.0 !wavenumber = wcm*1.0*w + monomer_E*1.0 !w = (w-1)*wcm*1.d0 + monomer_E*1.0 !w = w*wcm*1.d0 !w = w*monomer_E*1.0 !print *, w !wavenumber = w*wcm + monomer_E wavelength = 10000000.0/w !write(UNIT=4, FMT="(4(F15.7,2X))") Ab(x1), w !write(UNIT=4, FMT="(2(f15.10,2X))") w, Ab(x1) !write(UNIT=4, FMT="(4(F15.10,2X))") wavenumber, Ab(x1) !write(UNIT=4, FMT=105) weV, AbX(x1), AbY(x1), Abtot(x1) !here for eV! write(UNIT=4, FMT=105) weV, AbX(x1), AbY(x1), Abtot(x1) !here for wavelength! ! write(UNIT=4, FMT=105) w, AbX(x1), AbY(x1), Abtot(x1) !here for cm-1! 105 format(4e20.12) !print *, 'w and abs <',w,'/',Ab(x1),'>' end do close(UNIT=4) print *, 'kount1=',kount1,'' print *, 'kount2=',kount2,'' print *, 'kount3=',kount3,'' print *, 'kount4=',kount4,'' print *, 'kount5=',kount5,'' print *, 'Dimension of Hamiltonian=',B+1,'' print *, 'This is Eignvalue matrix' !---------------------------------------------------------CT Character-------------------------------------------------------------! print *, 'Here is the Vector' open(UNIT=6,FILE="LOOK5.txt",STATUS="OLD",ACTION="READWRITE") do x1=1,B wavelength = (10000000.0/Freq(x1)) write(UNIT=6, FMT=106) Wavelength, (COF1PE(x1)+COF2PE(x1)+COF3PE(x1)), (COFCTnn(x1)+COFCTnnv(x1)), COFDpAm(x1) 106 format(4e20.12) !print *, 'w and abs <',w,'/',Ab(x1),'>' end do close(UNIT=6) !-----------------------------------------------------------------------------------------------------------------------------------------! !open(UNIT=5,FILE="EGN1.txt",STATUS="OLD",ACTION="READWRITE") ! do i1=1,B ! write(UNIT=5, FMT="(2(F15.7,2X))") Eign(i1), F(i1) !write(*,322) Eign(i1), F(i1) !end do !322 format (6f10.1) ! close(UNIT=5) !-----------------------------------------------Dynamics-------------------------------------------------------------! !set the initial state: !--------------------------------------------------------------------------------------------------------------------------------------! end program Copolymer !----------------------------------------------------------------------------------------------! !-----------------------------------LL Diagonal Disorder--------------------------------------------! !--------------------------------------Disorder Subroutine------------------------------------------! Subroutine disorder_table() use common_variables implicit none !-------------------------------------! !integer :: Vx, Vy !Double Precision :: rand !real*8, external :: dlarnd !integer :: config, i, j, configmax, N !integer, parameter :: idist = 3 !Double Precision, dimension(:,:,:), allocatable :: disorder_elements !integer :: iseed(4)=(/47,3093,1041,77/) !-----------------------! call random_seed() do config=1,realization do Vx=1,N do Vy=1,N call random_number(rand1) call random_number(rand2) randR = dsqrt(-2*log(rand1)) theta = 2.0*PI*rand2 rand1N = randR*cos(theta) rand2N = randR*sin(theta) rand1Nnew = (rand1N*sigma) + randmean disorder_elements(config,Vx,Vy) = rand1Nnew end do end do end do !--------------------------! !to calculate the distribution! whole = 0 wholeEl = 0.0d0 do config=1,realization do Vx=1,N do Vy=1,N whole = whole + 1 wholeEl = wholeEl + disorder_elements(config,Vx,Vy) end do end do end do !------------------------------------! wholemean = wholeEl/whole !mean! wholedif = 0.0d0 do config=1,realization do Vx=1,N do Vy=1,N wholedif = wholedif + ((disorder_elements(config,Vx,Vy)-wholemean)**2) end do end do end do !-------------------------! !standarddeviation! standD = (1.0/(whole - 1))*wholedif standD = sqrt(standD) !---------------------------------! !Okay, we have the uniform distribution! open(UNIT=5,FILE="random3.txt",STATUS="NEW",ACTION="READWRITE") do config=1,realization do Vx=1,N do Vy=1,N !wholeP(config,Vx,Vy) = (1.0/sqrt(2.0*PI*(standD**2)))*EXP(-1.0*((disorder_elements(config,Vx,Vy)**2)/(2.0*standD**2))) wholeP(config,Vx,Vy) = (1.0/sqrt(2.0*PI*(sigma**2)))*EXP(-1.0*((disorder_elements(config,Vx,Vy)**2)/(2.0*sigma**2))) write(UNIT=5, FMT="(4(F15.7,2X))") disorder_elements(config,Vx,Vy), wholeP(config,Vx,Vy) end do end do end do close(UNIT=5) print *, wholemean, standD end Subroutine disorder_table !--------------------------------------------------------------------------------------------! !-----------------------------Auxilary routine for PRINT_MATRIX Subroutine----------------------! Subroutine PRINT_MATRIX(DESC, M, N, A, LDA) character :: DESC integer :: M, N, LDA Double Precision :: A(LDA, *) integer :: i, j Write(*,*) Write(*,*) DESC do i=1,M Write(*,9998) (A(i,j), j=1,N) end do 9998 format(11(:,1X, F6.2)) return end Subroutine PRINT_MATRIX !------------------------------------Coulomb Repulsion-----------------------------------------------------------! Function Repulsion(m,n,distance) implicit none integer :: m, n Double Precision :: distance, Repulsion, Rcoulomb Rcoulomb = iabs(m-n)*distance Repulsion = 8.3/Rcoulomb print *, 'Repulsion for <',m,'/',n,'> is ',repulsion,'' end Function Repulsion !---------------------------------------------Dopant-------------------------------------------------------------! Function Dopant(m,n,totaln,distance,danion) implicit none integer :: m, n, totaln, jm, jn Double Precision :: distance, Dopant, Rcoulombm, danion, Rcoulombn Double Precision :: dpm, dpn !notice, I have assumed that we have even number of chromophores! if (m.LE.(INT(totaln/2)).and.n.LE.(INT(totaln/2))) then dpm = ((INT(totaln/2)-m)*1.0*distance)+(distance/2.0) dpn = ((INT(totaln/2)-n)*1.0*distance)+(distance/2.0) Rcoulombm = sqrt((dpm*dpm)+(danion*danion)) Rcoulombn = sqrt((dpn*dpn)+(danion*danion)) Dopant = (-1.0*(8.3/Rcoulombm)) + (-1.0*(8.3/Rcoulombn)) else if (m.GT.(INT(totaln/2)).and.n.LE.(INT(totaln/2))) then jm = m-INT(totaln/2)-1 dpm = (jm*1.0*distance)+(distance/2.0) dpn = ((INT(totaln/2)-n)*1.0*distance)+(distance/2.0) Rcoulombm = sqrt((dpm*dpm)+(danion*danion)) Rcoulombn = sqrt((dpn*dpn)+(danion*danion)) Dopant = (-1.0*(8.3/Rcoulombm)) + (-1.0*(8.3/Rcoulombn)) else if (m.LE.(INT(totaln/2)).and.n.GT.(INT(totaln/2))) then jn = n-INT(totaln/2)-1 dpm = ((INT(totaln/2)-m)*1.0*distance)+(distance/2.0) dpn = (jn*1.0*distance)+(distance/2.0) Rcoulombm = sqrt((dpm*dpm)+(danion*danion)) Rcoulombn = sqrt((dpn*dpn)+(danion*danion)) Dopant = (-1.0*(8.3/Rcoulombm)) + (-1.0*(8.3/Rcoulombn)) else jm = m-INT(totaln/2)-1 jn = n-INT(totaln/2)-1 dpm = (jm*1.0*distance)+(distance/2.0) dpn = (jn*1.0*distance)+(distance/2.0) Rcoulombm = sqrt((dpm*dpm)+(danion*danion)) Rcoulombn = sqrt((dpn*dpn)+(danion*danion)) Dopant = (-1.0*(8.3/Rcoulombm)) + (-1.0*(8.3/Rcoulombn)) end if !print *, 'Dopant for <',m,'/',n,'> is ',Dopant,'' end Function Dopant !------------------------------------------------------CT Energy Function----------------------------------------! Function EnergyCT(m,n,wdnap,wdpan,Ntot,Sfactor,VMU) implicit none integer :: m, n, t, Ntot Double Precision :: EnergyCT, wdnap, wdpan Double Precision :: Sfactor, VMU !Double Precision :: nA, nD1, nD2 if ((MOD(m,2).eq.0).and.(iabs(m-n).EQ.1)) then EnergyCT = (Sfactor*wdnap) + VMU !print *, 'ECT <',m,'/',n,'> is ',EnergyCT,'' else if ((MOD(m,2).eq.0).and.((m-n).EQ.Ntot-1)) then EnergyCT = (Sfactor*wdnap) + VMU else if ((MOD(m+1,2).eq.0).and.(iabs(m-n).EQ.1)) then EnergyCT = (Sfactor*wdpan) + VMU !print *, 'ECT <',m,'/',n,'> is ',EnergyCT,'' else if ((MOD(m+1,2).eq.0).and.((m-n).EQ.1-Ntot)) then EnergyCT = (Sfactor*wdpan) + VMU end if end Function EnergyCT !------------------------------------------------------Te Function----------------------------------------! Function Tefunction(m,n,Teintra,Ntot) implicit none integer :: m, n, Ntot Double Precision :: Tefunction, Teintra !Double Precision :: nA, nD1, nD2 if ((MOD(m,2).eq.0).and.(iabs(n-m).EQ.1)) then Tefunction = Teintra !print *, 'Tef is <',m,'/',n,'> is ',Tefunction,'' else if ((MOD(m+1,2).eq.0).and.(iabs(n-m).EQ.1)) then Tefunction = Teintra !print *, 'Tef is <',m,'/',n,'> is ',Tefunction,'' else if ((MOD(m+1,2).eq.0).and.((m-n).EQ.1-Ntot)) then Tefunction = Teintra !print *, 'Tef is <',m,'/',n,'> is ',Tefunction,'' else if ((MOD(m,2).eq.0).and.((m-n).EQ.Ntot-1)) then Tefunction = Teintra !print *, 'Tef is <',m,'/',n,'> is ',Tefunction,'' end if end Function Tefunction !------------------------------------------------------Th Function----------------------------------------! Function Thfunction(m,n,Thintra,Ntot) implicit none integer :: m, n, Ntot Double Precision :: Thfunction, Thintra !Double Precision :: nA, nD1, nD2 if ((MOD(m,2).eq.0).and.(iabs(n-m).EQ.1)) then Thfunction = Thintra !print *, 'Thf is <',m,'/',n,'> is ',Thfunction,'' else if ((MOD(m+1,2).eq.0).and.(iabs(n-m).EQ.1)) then Thfunction = Thintra !print *, 'Thf is <',m,'/',n,'> is ',Thfunction,'' else if ((MOD(m+1,2).eq.0).and.((m-n).EQ.1-Ntot)) then Thfunction = Thintra !print *, 'Thf is <',m,'/',n,'> is ',Thfunction,'' else if ((MOD(m,2).eq.0).and.((m-n).EQ.Ntot-1)) then Thfunction = Thintra !print *, 'Thf is <',m,'/',n,'> is ',Thfunction,'' end if end Function Thfunction !----------------------------------------------Jfunction--------------------------------------------------------! Function Jfunction(m,n,JDA,JDD,JAA,Ntot) implicit none integer :: m, n, Ntot Double Precision :: Jfunction, JDA, JDD, JAA !Double Precision :: nA, nD1, nD2 !print *, ' JDA ',JDA,'' if ((MOD(m,2).eq.0).and.(iabs(n-m).EQ.1)) then Jfunction = JDA else if ((MOD(m+1,2).eq.0).and.(iabs(n-m).EQ.1)) then Jfunction = JDA else if ((MOD(m,2).eq.0).and.((n-m).EQ.1-Ntot)) then Jfunction = JDA else if ((MOD(m+1,2).eq.0).and.((n-m).EQ.Ntot-1)) then Jfunction = JDA else if ((MOD(m,2).eq.0).and.(iabs(n-m).EQ.2)) then Jfunction = JAA else if ((MOD(m+1,2).eq.0).and.(iabs(n-m).EQ.2)) then Jfunction = JDD !else if ((MOD(m,2).eq.0).and.((n-m).EQ.Ntot-2)) then !Jfunction = JAA !else if ((MOD(m+1,2).eq.0).and.((n-m).EQ.Ntot-2)) then !Jfunction = JDD end if end Function Jfunction !----------------------------------------Dipoles in different directions------------------------------------------! Function Dipoles(m,n,MuA,MuD,Angle,Alfadegree) implicit none integer :: m,n Double Precision :: Dipoles, MuA, MuD Double Precision :: Angle, Alfadegree Double Precision :: DegtoRad, MuxA, MuyA, MuxD, MuyD Real, parameter :: PI = 3.1415927 DegtoRad = Angle*PI/Alfadegree !Angle in radian! MuxA = MuA*sin(DegtoRad) MuyA = MuA*cos(DegtoRad) MuyD = MuD !Notice that the angle is expressed with respect to y-axis along which the MuD is poiting! MuxD = 0.0 !print *, ' MuA ',MuA,'' !print *, ' cosAng ',cos(DegtoRad),'' !print *, ' SinAng ',sin(DegtoRad),'' if ((m.EQ.1).and.(MOD(n,2).eq.0)) then !in y-direction! Dipoles = MuyA else if ((m.EQ.1).and.(MOD(n+1,2).eq.0)) then Dipoles = MuyD else if ((m.EQ.2).and.(MOD(n,2).eq.0)) then !9n x-direction! Dipoles = MuxA else if ((m.EQ.2).and.(MOD(n+1,2).eq.0)) then Dipoles = MuxD end if end Function Dipoles !-------------------------------------Frank-Condon Factor, general formula---------------------------------------! Function FC(m,n,lam) implicit none integer :: i, j, k, m, n !m: #quanta in ground state n: #quanta in excited state! Double Precision :: S, FC, C1, D1, Summ, C2, D2, C3, lam, DEN Double Precision, External :: Factorial S = lam*lam !lam is HR factor! C1 = dsqrt(Factorial(m)) D1 = dsqrt(Factorial(n)) !j = MIN(m,n) j = n Summ = 0.0 do i=1,(j+1) k = i - 1 if ((m-n+k).lt.0) cycle C2 = Factorial(k) D2 = Factorial((n-k)) C3 = Factorial((m-n+k)) DEN = C2*D2*C3 Summ = Summ + ((((-1)**(m-n+k))*(S**((m-n+(2*k))/2.0)))/DEN) end do FC = C1*D1*dexp(-1.0*(S/2.0))*Summ !print *, 'FC factor for <',m,'/',n,'> is ',FC,'' !print *, 'lamda ',lam,'' end Function FC !-------------------------------------------------------------------------------------------------------------------! !-------------------------------------Factorial--------------------------------------------------------------------! Function Factorial(x) implicit none integer :: i, x Double Precision :: Factorial Factorial = 1.0 do i=1,x Factorial = Factorial*(1.0*i) end do !print *, 'Factorial of ',x,' is ',Factorial,'' end Function Factorial !-------------------------------------------------------------------------------------------------------------!