*******************************************************************************

	program pesgen_thermal 

*******************************************************************************
c
c  Program to extract binding energies and grand canonical corrections
c  from a series of files.
c
c	To compile non-graphics programs use : 
c             f77 $(OPT) pesgen.f -o pesgen 
c
c
	character*72    atomform
	character*80    aform 
	character*80    record, rec
	character*24    mulform
	character*80    freqform
        character*40	fname,gname
        character*40	hname
        character*4     root,root2
        character*1     tchar1,tempchar1
        character*2     tchar2,tempchar2
	real*8          rdist(100)
	real*8          kcal(100)
	real*8          ev(100)
	real*8          au(100)
	real*8          therm(100),th
	integer  	iatom,im,i,j,ierror
	integer  	ii,il,ir
	integer		iunit 
	integer		itherm
c
c 
         itherm=1 
c
	 iunit=3
	 print *, 'Enter first outmol file rootname: '
	 read(5,'(a4)') root2
c
	 print *,'Enter the number of outmol files: '
         read(5,*) kfiles 
c
	 print *,'Enter starting distance: '
         read(5,*) rstart 
	 print *,'Enter increment: '
         read(5,*) rinc
         ii = 1
         do 700 k=1,kfiles
          kk=k
          call int2char(kk,tchar1,tchar2)
            if (kk.le.9) then
               fname = root2//tchar1//'.outmol'
            else
               fname = root2//tchar2//'.outmol'
            endif
       open(unit=3,name=fname,status='old',err=300)
c
        print *,'Reading file ',fname,'...'
        do 200 j=1,30000
        read(3,'(a80)',end=400,err=999) record
c
c find the grand canonical ensemble correction
c
        if(itherm.eq.1) then
        if(record(1:12).eq.'Fermi Energy') then
          read (record,42,err=445) therm(ii)
334     continue
          if(therm(ii).gt.0.0) then
            therm(ii)=-therm(ii)
            Print *,'Something may be wrong with the format'
         endif
        goto 554
445     continue
        read (record,43,err=997) therm(ii)
        goto 334
554     continue
        endif
        endif
        if(record(1:14).eq.'binding energy') then
c        print *,record
        read (record,30,err=999) au(ii),ev(ii),kcal(ii)
c       print *,kcal(ii)
        goto 500
        endif
200     continue
c
         goto 400
c
300     continue
         print *,'Could not find file ',fname
         goto 600
400     continue 
         print *,'No data found in file ',fname
         goto 550
500     continue
         rdist(ii)=rstart+(k-1)*rinc
         ii=ii+1
550     continue
         close(3)
600     continue 
700     continue
c
	 print *, 'Enter second outmol file rootname: '
	 read(5,'(a4)') root
c
	 print *,'Enter the number of outmol files: '
         read(5,*) kfiles 
c
	 print *,'Enter starting distance: '
         read(5,*) rstart 
	 print *,'Enter increment: '
         read(5,*) rinc
         do 2700 k=1,kfiles
          call int2char(k,tchar1,tchar2)
            if (k.le.9) then
               fname = root//tchar1//'.outmol'
            else
               fname = root//tchar2//'.outmol'
            endif
       open(unit=3,name=fname,status='old',err=2300)
c
c        print *,'Reading file ',fname,'...'
        do 900 j=1,30000
        read(3,'(a80)',end=2400,err=999) record
        if(itherm.eq.1) then
        if(record(1:12).eq.'Fermi Energy') then
          read (record,42,err=2445) therm(ii)
2334     continue
          if(therm(ii).gt.0.0) then
            therm(ii)=-therm(ii)
            Print *,'Something may be wrong with the format'
         endif
        goto 2554
2445     continue
c        print *, record
        read (record,43,err=997) therm(ii)
        goto 2334
2554     continue
         endif
        endif
        if(record(1:14).eq.'binding energy') then
c        print *,record
        read (record,30,err=999) au(ii),ev(ii),kcal(ii)
c       print *,kcal(ii)
        goto 2500
        endif
900     continue
c
         goto 2400
c
2300     continue
         print *,'Could not find file ',fname
         goto 2600
2400     continue 
         print *,'No data found in file ',fname
         goto 2550
2500     continue
         rdist(ii)=rstart+(k-1)*rinc
         ii=ii+1
2550     continue
         close(3)
2600     continue 
2700     continue
c
         gname=root2//'.output'
         junit=7
         open(unit=7,file=gname,form='formatted',err=999)
         do 100 k=1,ii-1
          th=therm(k)*2625.55
          energy=kcal(k)*4.184-therm(k)*2625.55
          write (junit,25) rdist(k),kcal(k)*4.184,th
100     continue
        close(7)
         gname=root2//'.kj'
         junit=7
         open(unit=7,file=gname,form='formatted',err=999)
         energyref=kcal(ii-1)*4.184-therm(ii-1)*2625.55
         do 170 k=1,ii-1
          if(itherm.eq.1) then
            energy=kcal(k)*4.184-therm(k)*2625.55
            energy=energy-energyref
            write (junit,20) rdist(k),energy
          else
            write (junit,20) rdist(k),kcal(k)*4.184
          endif
170     continue
        close(7)
c
20      format(f10.5,' , ',f14.5)
25      format(f10.5,' , ',f14.5,' , ',f14.5)
30      format(14x,f11.7,2x,f11.5,2x,f11.3,9x)
42      format(48x,f10.7,22x)
43      format(46x,f10.7,24x)
        goto 3001
997     continue
        print *,'Error in Fermi Energy line'
        goto 3002
999      continue
        print *,'Error in file'
        goto 3002
3001    continue
        print *,'PROGRAM COMPLETE'
3002    continue
	end

***************************************************************************

      Subroutine int2char(i,tempchar1,tempchar2)

***************************************************************************
!          To convert the integer data to character data

      Implicit None

      Integer i
      Character*1 tempchar1
      Character*2 tempchar2
      Character*4 record
 
      if(i.lt.10) then
        Write (record,'(i1)') i
        Read (record,'(a1)') tempchar1
      else
        Write (record,'(i2)') i
        Read (record,'(a2)') tempchar2
      endif

      Return
      End