program pothist
c     histogram from potential energy distribution
      implicit none
      integer*4 nc,ic,it,is,i,j,im,iw,np,ncoord
      real*8 w,wmin,wmax,dw,an
      character*80 s80
      logical lex
      real*8,allocatable::h(:,:),p(:)
      integer*4,allocatable::s(:,:)
      character*1 norm

      ncoord=15

c     1 - strecth
c     2 - bend
c     3 - torsion
c     4 - OPLA

c     5 - IM, stretch
c     6 - IM, bend, 4 types
c     7 - IM, torsion

c     8 - TX, helix
c     9 - TY and TZ, helix
c    10 - RX, helix
c    11 - RY and RZ, helix

c    12 - H-stretch
c    13 - H-breath
c    14 - H-def
c    15 - H-unwinding

      wmin=0.0d0
      wmax=4000.0d0
      np=400
      dw=(wmax-wmin)/dble(np)
      allocate(h(ncoord,np),s(ncoord,2))

      inquire(file='HIST.TXT',exist=lex)
      if(lex)then
       open(9,file='HIST.TXT')
       do 8 i=1,np
       w=wmin+(i-1)*dw
8      read(9,*)h(1,i),(h(j,i),j=1,ncoord)
       close(9)
      else
       do 5 j=1,ncoord
       do 5 i=1,np
5      h(j,i)=0.0d0
      endif

      norm='n'
      if(iargc().gt.0)call getarg(1,norm)
      if(norm.eq.'y')norm='Y'
      if(norm.eq.'Y')then
c      just normalize histograms and quit:
       do 9 i=1,ncoord
       an=0.0d0
       do 10 j=1,np
10     an=an+h(i,j)
       if(an.gt.1.0d-5)then
        do 11 j=1,np
11      h(i,j)=h(i,j)/an/dw
       endif
9      continue
       write(6,*)'Histograms have been normalized'
       goto 77
      endif

      do 2 i=1,ncoord
      do 2 j=1,2
2     s(i,j)=0
      open(9,file='FILE.UMA')
      read(9,*)
      read(9,*)
      read(9,*)nc,nc
      ic=0
1     read(9,90,end=99,err=99)s80
90    format(a80)
      if(s80(56:62).eq.'STRETCH')then
       ic=ic+1
       if(s(1,1).eq.0)s(1,1)=ic
       s(1,2)=ic
      endif
      if(s80(56:59).eq.'BEND')then
       ic=ic+1
       if(s(2,1).eq.0)s(2,1)=ic
       s(2,2)=ic
      endif
      if(s80(56:62).eq.'TORSION')then
       ic=ic+1
       if(s(3,1).eq.0)s(3,1)=ic
       s(3,2)=ic
      endif
      if(s80(56:59).eq.'OPLA')then
       ic=ic+1
       if(s(4,1).eq.0)s(4,1)=ic
       s(4,2)=ic
      endif
      if(s80(56:62).eq.'INTERMO')then
       ic=ic+1
       read(s80(1:6),*)it
       if(it.eq.8)is=5
       if(it.eq.9)is=6
       if(it.eq.10)is=6
       if(it.eq.11)is=6
       if(it.eq.12)is=6
       if(it.eq.13)is=7
       if(s(is,1).eq.0)s(is,1)=ic
       s(is,2)=ic
      endif

      if(s80(56:64).eq.'H-TX')then
       ic=ic+1
       if(s(8,1).eq.0)s(8,1)=ic
       s(8,2)=ic
      endif
      if(s80(56:64).eq.'H-TY')then
       ic=ic+1
       if(s(9,1).eq.0)s(9,1)=ic
       s(9,2)=ic
      endif
      if(s80(56:64).eq.'H-TZ')then
       ic=ic+1
       if(s(9,1).eq.0)s(9,1)=ic
       s(9,2)=ic
      endif

      if(s80(56:64).eq.'H-RX')then
       ic=ic+1
       if(s(10,1).eq.0)s(10,1)=ic
       s(10,2)=ic
      endif
      if(s80(56:64).eq.'H-RY')then
       ic=ic+1
       if(s(11,1).eq.0)s(11,1)=ic
       s(11,2)=ic
      endif
      if(s80(56:64).eq.'H-RZ')then
       ic=ic+1
       if(s(11,1).eq.0)s(11,1)=ic
       s(11,2)=ic
      endif

      if(s80(56:64).eq.'H-stretch')then
       ic=ic+1
       if(s(12,1).eq.0)s(12,1)=ic
       s(12,2)=ic
      endif
      if(s80(56:64).eq.'H-breath.')then
       ic=ic+1
       if(s(13,1).eq.0)s(13,1)=ic
       s(13,2)=ic
      endif
      if(s80(56:61).eq.'H-def.')then
       ic=ic+1
       if(s(14,1).eq.0)s(14,1)=ic
       s(14,2)=ic
      endif
      if(s80(56:63).eq.'H-unwind')then
       ic=ic+1
       if(s(15,1).eq.0)s(15,1)=ic
       s(15,2)=ic
      endif
      goto 1
99    close(9)
      do 3 i=1,ncoord
3     write(6,600)i,s(i,1),s(i,2)
600   format(i4,':',i4,' - ',i4)

      if(ic.ne.nc)then
       write(6,*)ic,nc
       write(6,*)'Some coordinates missed'
       stop
      endif

      allocate(p(nc))
      open(9,file='POT_all.LST')
      do 4 im=1,nc
      read(9,*,end=88,err=88)w,w,(p(i),i=1,nc)
      if(w.lt.10.0d0)goto 88
      iw=int((w-wmin)/dw)
      if(iw.ge.1.and.iw.le.np)then
       do 6 i=1,nc
       do 6 j=1,ncoord
6      if(i.ge.s(j,1).and.i.le.s(j,2))h(j,iw)=h(j,iw)+p(i)
      endif
 4    continue
 88   close(9)

77    open(9,file='HIST.TXT')
      do 7 i=1,np
      w=wmin+(i-1)*dw
      write(9,70)w
70    format(f12.2,$)
      do 71 j=1,ncoord
71    write(9,73)h(j,i)
73    format(g12.4,$)
7     write(9,*)
      close(9)

      end