      program decoy_set_analysis_for_orient_potential
!     calculating decoys from potential
      use protein_information
      use energy_function
      implicit none
      integer,parameter :: nmax=(1000)

      character*255 dir_decoy,dir_energy,list_decoy,dir_result
      integer num,eof,rank,num1,num_decoy
      real*8 energy,rmsd,native
      character*255 pdblist,filename,pdbname,subset
      character*255 a1,a2,a3,a4,a5,a6
      type(residue_type),dimension(nmax)::residues,res_native,res_decoy
     c     ,res_native0
      integer i

      call read_input(dir_energy,dir_decoy,list_decoy,dir_result)

      call init(dir_energy)

      open(unit=2,file=trim(dir_decoy)//trim(list_decoy),status='old')
      do
         read(2,*,IOSTAT=eof) subset
         if(eof.lt.0)exit
         print *,trim(subset)," subset:"

         open(unit=1,file=trim(dir_decoy)//trim(subset)//"/list"
     c        ,status='old')
         do
            read(1,*,IOSTAT=eof) pdbname
            if(eof.lt.0)exit

            open(unit=20,file=trim(dir_decoy)//trim(subset)//"/"//
     c           trim(pdbname)//"/list",status='old')
!            open(unit=21,file=trim(dir_decoy)//trim(subset)//"/"//
!     c           trim(pdbname)//"/rmsds",status='old')

            open(unit=31,file=trim(dir_result)//trim(subset)//"_"
     c           //trim(pdbname),status='unknown')

            rank=1
            pdblist=trim(pdbname)//".pdb"
            rmsd=0.

            filename=trim(dir_decoy)//trim(subset)//"/"
     c           //trim(pdbname)//"/"//pdblist

            call import_pdb(nmax,filename,num,res_native)
            
            print *,"Analyzing on ",trim(pdbname),"(",num," residues)"

!            print "(5000I3)",(res_native(i)%res_type,i=1,num)

            read(20,*)pdblist
            do while(trim(pdblist).eq.trim(pdbname)//".pdb")
               read(20,*)pdblist
            enddo
            filename=trim(dir_decoy)//trim(subset)//"/"//
     c           trim(pdbname)//"/"//pdblist

            call import_pdb(nmax,filename,num_decoy,res_decoy)
            
!            print "(A4,2I6,F10.5)",trim(pdbname),num,
!     c           num_decoy,num_decoy*1./num

            call native_vs_decoy(num,num_decoy,res_native,res_decoy,
     c           res_native0)

!            call export_pdb(num_decoy,res_native0,"test.pdb")
            
!            open(unit=37,file=trim(pdbname)//'_e.dat',
!     c           status='unknown')

            call ca_energy(num_decoy,res_native0,energy)

!            close(37)
!            print *,energy
!            stop

            write(31,"(2F15.4)")rmsd,energy

            native=energy
            eof=0
            rewind(20)
!            cycle
            do
               read(20,*,IOSTAT=eof) pdblist
               if(eof.lt.0)exit
!               read(21,*,IOSTAT=eof) a1,a2,a3,a4,a5,a6,rmsd
!               if(eof.lt.0)exit

               if(trim(pdblist).ne.trim(pdbname)//".pdb")
     c              then
                  filename=trim(dir_decoy)//trim(subset)//"/"//
     c                 trim(pdbname)//"/"//pdblist

                  call import_pdb(nmax,filename,num1,residues)

                  if(num_decoy.ne.num1)then
                     print *,"decoy protein size is not consistant"
                     stop
                  endif

!                  open(unit=37,file=trim(pdblist)//'_e.dat',
!     c                 status='unknown')

                  call ca_energy(num_decoy,residues,energy)

!                  close(37)

                  call protein_rmsd(num_decoy,res_native0,residues,
     c                 rmsd,3,1)

                  write(31,"(2F15.4)")rmsd,energy
                  if(energy.lt.native)rank=rank+1

               endif
            enddo
            close(20)
!            close(21)
            close(31)
            print *,"Rank of ",trim(pdbname)," native: ",rank
         enddo
         close(1)
      enddo
      close(2)
      end
c
      subroutine read_input(dir_energy,dir_decoy,list_decoy,dir_result)
      implicit none 
      character*255 dir_energy,dir_decoy,list_decoy,dir_result,a1
      logical ex

      inquire(file='input.dat',exist=ex)
      if(.not.ex)then
         print *,"Input file not exists!"
         call system("echo energy = ./energy/ >> input.dat")
         call system("echo decoy = ./ >> input.dat")
         call system("echo list = list>> input.dat")
         call system("echo result = ./ >> input.dat")

         print *,"an example input file is generated!"
         call system("cat input.dat")
         print *,"press any key to continue...."
         read(*,*)
      endif

      open(unit=1,file='input.dat',status='old')
      read(1,"(A9,A255)")a1,dir_energy
      read(1,"(A8,A255)")a1,dir_decoy
      read(1,"(A7,A255)")a1,list_decoy
      read(1,"(A9,A255)")a1,dir_result
      close(1)

      return
      end subroutine read_input
c
      subroutine native_vs_decoy(num,num_decoy,res_native,res_decoy,
     c     res_native0)
      use protein_information
      implicit none
      integer num,num_decoy
      type(residue_type)res_native(num),res_decoy(num_decoy),
     c     res_native0(num_decoy)
      integer i,j,trail,head
      logical match

!      print *,(res_native(i)%res_type,i=1,num)
!      print *,""
!      print *,(res_decoy(i)%res_type,i=1,num_decoy)

      trail=num-num_decoy+1

      match=.false.
      do i=1,trail         
         do j=1,num_decoy
            if(res_native(j+i-1)%res_type.ne.res_decoy(j)%res_type)exit
            if(j.eq.num_decoy)match=.true.
         enddo         
         if(match)then
            head=i
            exit
         endif
      enddo

      if(.not.match)then
         print *,"error, no match found"
         stop
      else
         do i=1,num_decoy
            res_native0(i)=res_native(i+head-1)
         enddo
      endif

      return
      end subroutine native_vs_decoy
