Code:
program check
implicit none
! Setting up parameters and necessary conversion factors
real,parameter::h=4.135667516e-15; ! planck's constant in eVs
real,parameter::c=3.0e10; ! velocity of light in cm/s
real,parameter::kb=8.6173324e-5; ! boltzmann's constant in eV/K
integer, parameter:: dp = selected_real_kind(15, 307)
real(dp),allocatable::zper(:) ! zero point energy array
real(dp),allocatable::zpcr(:) ! zero point corrected energy
real(dp),allocatable::qvibr(:)
real(dp)::zpcrtot=0.0 ! total zero point corrected energy
integer::i,j
integer::r
character(len=20)::filename
integer::ierror,alloc_err,nlines,nless,array_t
real(dp)::scf,vibmulti,freqtot
real(dp),allocatable::freq(:),vib(:)
real(dp)::T(29) =(/(array_t,array_t=273,973,25)/)
write(*,*)"How many reactants do you have?"
read(*,*)r
!allocate proper dimension to all zero point arrays
allocate(zper(r))
allocate(zpcr(r))
allocate(qvibr(r))
! Loop through reactant files
do j= 1,29 !temp loop
write(*,*)"This is for temperature: ",T(j)
do i=1,r ! file loop
write(filename,100)i
100 format('reac-',I1) ! get filenames
write(*,*)"This is for reac-",i
freqtot=0 ! initializing total frequency value
vibmulti=1 !
! Count no of lines
nlines=0
OPEN(unit=10,file=filename,status='old',action='read',&
iostat=ierror)
do ! count loop
read(10,*,iostat=ierror)
if (ierror /= 0) exit
nlines=nlines+1
end do !end count loop
CLOSE(unit=10)
! Count no of frequency elements
nless=nlines-1
!allocate a dimension to frquency array
allocate(freq(nless))
allocate(vib(nless))
OPEN(unit=10,file=filename,status='old',action='read',&
iostat=ierror)
read(10,*,iostat=ierror) scf
do ! freq loop
read(10,*,iostat=ierror) freq(nless)
if (ierror /= 0) exit
if (freq(nless) < 100) then
freq(nless) = 100
end if
freqtot=freqtot+freq(nless)
vib(nless)= 1/(1-exp(-h*c*freq(nless)/kb/T(j)))
vibmulti=vib(nless)*vibmulti
end do !end freq loop
deallocate (freq,stat=alloc_err)
deallocate (vib,stat=alloc_err)
CLOSE(10)
zper(i)= 0.5*h*c*freqtot ! zero point energy
zpcr(i)= scf + zper(i) ! zero point correction
zpcrtot= zpcrtot+zpcr(i) ! total sum of zero point corrected energies of reactants
write(*,*)"ZPER", zper(i)
write(*,*)"ZPCR", zpcr(i)
write(*,*)"ZPCRTOT", zpcrtot
write(*,*)"Vib.P.F for reac",i,vibmulti
qvibr(i)=vibmulti
write(*,*)"This is test qvib: ", qvibr(i)
end do ! end file loop
write(*,*) "Try product: ",PRODUCT(qvibr)
deallocate (qvibr,stat=alloc_err)
deallocate (zper,stat=alloc_err)
deallocate (zpcr,stat=alloc_err)
end do ! end temp loop
end program check