      program decoy_set_analysis_for_orient_potential_hr
!     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,nativefile,filelabel
      character*255 a1,a2,a3,a4,a5,a6
      type(residue_type),dimension(nmax)::residues,res_native
      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)//"/Test150.txt"
     c        ,status='old')
         do
            read(1,*,IOSTAT=eof) pdbname
            if(eof.lt.0)exit

!            pdbname="1g10A"

            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)//".native.pdb"
            rmsd=0.

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

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

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

            call ca_energy(num,res_native,energy)

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

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

            native=energy
            eof=0            
!            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

!               pdblist="1g10A.305.coord.pdb"

               filelabel=trim(dir_decoy)//trim(subset)//"/"//
     c              trim(pdbname)//"/"//trim(pdbname)//".llabel.pdb"

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

               call import_hr(num,filelabel,filename,residues)
!               call import_pdb(nmax,filename,num1,residues)

               call ca_energy(num,residues,energy)

               call protein_rmsd(num,res_native,residues,
     c              rmsd,3,1)

!               print*,rmsd,energy
!               call export_pdb(num,residues,trim(dir_result)//"out.pdb")

!               stop

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

            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 import_hr(num,filelabel,filename,residues)
      use protein_information
      use common_sense
      use constant
      implicit none
      integer num
      type(residue_type) residues(num)
      integer atomindex
      integer i,j,k,eof,filter,nindex,preindex
      real*8 x0,y0,z0
      character filename*255,head*13,type*3,temp,filelabel*255
      character gap*2,res*3,gap2*1,chid*1,other*4,last*26
      character precode*1,alt_now*1,alt*1,icode*1

      open(unit=15,file=trim(filelabel),status='old')
      open(unit=16,file=trim(filename),status='old')
      read(16,*)temp
      i=0
      j=1
      eof=0
      do
         read(15,2001,IOSTAT=eof) head,type,gap,res,
     c        gap2,chid,nindex
         if(eof.lt.0)then
            exit
         endif
         if(head(1:6).eq."ENDMDL")exit
         if((head(1:4).eq."ATOM").and.(gap.ne."B"))then
            if((i.eq.0).and.(res(1:3).eq.'ACE'))cycle
            if(type(1:1).eq.'H')cycle

            if((j.eq.1).or.
     c           ((preindex.ne.nindex).or.(icode.ne.precode).or.
     c           (residues(i)%res_type.ne.index(res_list,res)/4+1)))then
!            if((j.eq.1).or.(preindex.ne.nindex))then
               i=i+1
               residues(i)%res_type=index(res_list,res)/4+1
               preindex=nindex
               precode=icode
               alt_now=" "
               do k=1,14
                  residues(i)%r(1,k)=large
               enddo
            endif

            if((alt_now.eq." ").and.(alt.ne." "))then
               alt_now=alt
            endif

            if((i.eq.1).and.(type.eq."NT"))type="N"

            select case(type)
            case ('AE1 ')
               type='OE1 '
            case ('AE2 ')
               type='NE2 '
            case ('AD1 ')
               type='OD1 '
            case ('AD2 ')
               type='ND2 '
            end select

            atomindex=index(atom_list(residues(i)%res_type),type)
            if(atomindex.ne.0)then
               read(16,*)x0,y0,z0
               residues(i)%r(1,atomindex/3+1)=x0
               residues(i)%r(2,atomindex/3+1)=y0
               residues(i)%r(3,atomindex/3+1)=z0
               j=j+1
            elseif((type.ne."OXT").and.(type(1:1).ne."H").and.
     c              (type(1:1).ne.'Q').and.(type(1:1).ne.'D').and.
     c              (type.ne."OT"))then
               print *,"unknown atom type in residue ",res,' : ',type
            endif
         endif         
      enddo
!     print *,natom,nres
 2001 format(A13,A3,A1,A3,A1,A1,I4)
      close(15)
      close(16)
      return
      end subroutine import_hr
