      program Doppioni

      implicit none
      character*7 pdbfile,nomefile
      character*11 file_inp
      integer NUM,NUM_STEP,base1,base2,base3,base4,errore_inp,errore_pdb
      integer nlinemax,ndmax,residuo1(2000),residuo2(2000)
      integer residuo3(2000),residuo4(2000)
      integer ncount,nca,i,j,k,nres,ii,jj,useless,molecul
      integer t,differenza,numero_residui_ATOM,numero_residui_HETATM
      integer numero_residui,differenza_ATOM,differenza_HETATM
      integer incremento_ATOM,incremento_HETATM,appoggio
      character*4 Form
      character*30 inutile
      character*16 file_corretto
      character*5 char_base1,char_base2,char_base3,char_base4
      logical condizione

      parameter(nlinemax=500000)
      parameter(ndmax=1000)
     
      real*8 tau(3),pos(3,nlinemax),occupancy,pos_useless(3,nlinemax)
      real*8 X_CM,Y_CM,Z_CM,Modulo, X,Y,Z,Mod,rgyr
      real*8 pos_CM(3,nlinemax)
c      real*8 X(ndmax),Y(ndmax),Z(ndmax),Mod(ndmax)
      integer num_res(nlinemax)
      character*6 atom
      character*3 resname
      character*4 atname,anam(nlinemax),renam(nlinemax) !HO CAMBIATO PERCHE' PRIMA ERANO CARATTERI DI 3!!!!!!!!!!!!!!!
      character*4 nam1,nam2,re1,re2,nam3
      character*1 chain,ch(nlinemax)
      character*6 string
      character*30 filename

      real*8 rmin,rmax,rd,r(ndmax),rr
      real*8 amax,amin,ad,a(ndmax),an,rr1,rr2,rr3,coss
      integer nd,na,caca(ndmax),cap(ndmax),pp(ndmax)  !nd per la distanza; na per l'angolo e per il diedro. Chiaramente usando na per due cose
      integer aa(ndmax),gg(ndmax),cc(ndmax),uu(ndmax) !diverse mi dovro ricordare di reinizializzare certe variabili collegate a na (vedi riga 476)
      integer gc(ndmax),ga(ndmax),gu(ndmax)
      integer ca(ndmax),cu(ndmax),au(ndmax)
      integer nn(ndmax),mm(ndmax),np(ndmax),nm(ndmax),pm(ndmax)


      real*8 r1(3),r2(3),r3(3),n1(3),n2(3),cosalfa,ss,scal(3)
      real*8 alfa(nlinemax)

c---------------------------------------------------------------------
c     set the the range of r
c---------------------------------------------------------------------

      rmin=0.
      rmax=30.
      nd=100

      rd=(rmax-rmin)/nd !DIMENSIONE DEL BEAD


c---------------------------------------------------------------------
c     initialize the vector caca(), used as basic vector
c     for the distance distributions
c---------------------------------------------------------------------
      do i=1,nd
         r(i)=i*rd+rmin
         caca(i)=0
      enddo



      do i=1,2000
         residuo1(i)=-1000000
         residuo2(i)=-1000000
         residuo3(i)=-1000000
         residuo4(i)=-1000000

      enddo

      condizione=.false.

      open(10,file='percentuali',status='old',err=1) !LA LABEL 1 E' ALLA FINE

c      do i=1,100000
 5       read(10,110,end=999) pdbfile,NUM !RICORDA CHE DEVI METTERE NELLA CARTELLA I FILE PDB
         if(errore_inp.gt.0) then
            print*,'Manca il file ',file_inp
         else
            continue
         endif

         if(errore_pdb.gt.0) then
            print*,'Manca il file ',pdbfile
         else
            continue
         endif

         NUM_STEP=NUM-1
         numero_residui=0
         numero_residui_HETATM=0
         numero_residui_ATOM=0
         incremento_ATOM=0
         incremento_HETATM=0
         if(pdbfile(7:7).eq.'.') then
          
            file_inp=pdbfile // 'inp'
            nomefile=pdbfile(1:6)
         else
            file_inp=pdbfile // '.inp'
            nomefile=pdbfile
         endif
         print*,nomefile,NUM
         print*,file_inp
         open(15,file=file_inp,status='old',err=5,iostat=errore_inp) !ANCHE I FILE .INP VANNO NELLA CARTELLA
         open(12,file=nomefile,status='old',err=5,iostat=errore_pdb)
         do i=1,5
            read(15,120) inutile
         enddo


         read(15,200,end=4) char_base1,char_base2
        
         if((char_base1(5:5).eq.'_').and.(char_base2(5:5).eq.'_')) then
           
           if((char_base1(1:3).eq.'...').and.(char_base2(1:3).eq.'...')) 
     $           then
              backspace 15
              read(15,300,end=4) base1,base2
           else
    
          
           if((char_base1(1:3).eq.'...').and.(char_base2(1:2).eq.'..')) 
     $          then
              backspace 15
              read(15,400,end=4) base1,base2
           else
          
           
           if((char_base1(1:3).eq.'...').and.(char_base2(1:1).eq.'.')) 
     $          then
              backspace 15
              read(15,500,end=4) base1,base2
           else

                       
           if((char_base1(1:2).eq.'..').and.(char_base2(1:3).eq.'...')) 
     $          then
              backspace 15
              read(15,550,end=4) base1,base2
           else 
          
          
          if((char_base1(1:2).eq.'..').and.(char_base2(1:2).eq.'..')) 
     $         then
             backspace 15
             read(15,600,end=4) base1,base2
          else
            
          
          if((char_base1(1:2).eq.'..').and.(char_base2(1:1).eq.'.')) 
     $         then
             backspace 15
             read(15,700,end=4) base1,base2
          else

          if((char_base1(1:1).eq.'.').and.(char_base2(1:3).eq.'...')) 
     $            then
             backspace 15
             read(15,730,end=4) base1,base2
          else

          if((char_base1(1:1).eq.'.').and.(char_base2(1:2).eq.'..')) 
     $            then
             backspace 15
             read(15,760,end=4) base1,base2
          else
            
          
          if((char_base1(1:1).eq.'.').and.(char_base2(1:1).eq.'.')) 
     $         then
             backspace 15
             read(15,800,end=4) base1,base2
          else
             continue
          endif
       endif
      endif
      endif
      endif
      endif
      endif
      endif
      endif

         else
            continue
         endif



c 6000    read(15,130,err=6000) base1,base2
       
         do j=1,NUM_STEP
            read(10,120,end=4) Form
          
c 7000       read(15,130,err=7000,end=4) base3,base4 

          
        read(15,200,end=4) char_base3,char_base4
        
          if((char_base3(5:5).eq.'_').and.(char_base4(5:5).eq.'_')) then

           if((char_base3(1:3).eq.'...').and.(char_base4(1:3).eq.'...')) 
     $           then
              backspace 15
              read(15,300,end=4) base3,base4
           else
    
          
           if((char_base3(1:3).eq.'...').and.(char_base4(1:2).eq.'..')) 
     $          then
              backspace 15
              read(15,400,end=4) base3,base4
           else
          
           
           if((char_base3(1:3).eq.'...').and.(char_base4(1:1).eq.'.')) 
     $          then
              backspace 15
              read(15,500,end=4) base3,base4
           else

           if((char_base3(1:2).eq.'..').and.(char_base4(1:3).eq.'...')) 
     $             then
              backspace 15
              read(15,550,end=4) base3,base4
           else 
          
          
          if((char_base3(1:2).eq.'..').and.(char_base4(1:2).eq.'..')) 
     $         then
             backspace 15
             read(15,600,end=4) base3,base4
          else
            
          
          if((char_base3(1:2).eq.'..').and.(char_base4(1:1).eq.'.')) 
     $         then
             backspace 15
             read(15,700,end=4) base3,base4
          else
            
          if((char_base3(1:1).eq.'.').and.(char_base4(1:3).eq.'...')) 
     $            then
             backspace 15
             read(15,730,end=4) base3,base4
          else

          if((char_base3(1:1).eq.'.').and.(char_base4(1:2).eq.'..')) 
     $            then
             backspace 15
             read(15,760,end=4) base3,base4
          else

          if((char_base3(1:1).eq.'.').and.(char_base4(1:1).eq.'.')) 
     $         then
             backspace 15
             read(15,800,end=4) base3,base4
          else
            continue
         endif
      endif
      endif
      endif
      endif
      endif
      endif
      endif
      endif

      else
         continue
      endif

          
            if(Form.eq.'  A ') then
             
               residuo1(j)=base1
               residuo2(j)=base2
               residuo3(j)=base3 !POTREI AVERE PIU' DI UN ' A '
               residuo4(j)=base4
            else
               continue
            endif
            base1=base3
            base2=base4
               
         enddo

       
 4       do j=1,NUM_STEP
c            print*,'*********'
c            print*,residuo1(j),residuo2(j)
c            print*,residuo3(j),residuo4(j)
         enddo
        
c----------------------------
c una volta memorizzati i residui che mi servono per il tipo A
c procedo a calcolare le distanze
c Calclo tutte le distanze e poi seleziono solamente quelle che mi servono
c---------------------------------------

         nca=0

         do i=1,nlinemax        !QUESTO CICLO SERVE PER L'ISTRUZIONE READ SEGUENTE AFFINCHE IL FILE VENGA LETTO PER INTERO,
                      !SI INTENDE DALLA RIGA 1 ALLA RIGA LINEMAX.iNFATTI L'ISTRUZIONE DI READ SUCCESSIVA LEGGE UN
                      !PARTICOLARE FORMATO CHE E' IL 1000 QUINDI UNA VOLTA FATTA LA LETTURA IN QUESTO FORMATO LA 
                      !LETTURA SI ARRESTA ED HO BISOGNO DI UN CICLO PER FARLA PROSEGUIRE ALLA RIGA SUCCESSIVA DELLO STESSO FILE

 6       read(12,1000,err=6,end=11) string,ncount,atname,resname,chain,
     $        nres,(tau(k),k=1,3),occupancy
         
c         print*,string,ncount,atname,resname,chain,
c     $        nres,(tau(k),k=1,3),occupancy

c         if((atname.eq.'CA ').or.(atname.eq.'CA1').or.
c     $        (atname.eq.'CA2').or.(atname.eq.'CA3')
c     $        .or.atname.eq.'P  ') then
         
         if((nres.ne.0).and.((string(1:4).eq.'ATOM').or.(string.eq.
     $        'HETATM')).and.((resname.eq.'  A').or.(resname.eq.'  C')
     $        .or.(resname.eq.'  G').or.(resname.eq.'  T').or.(resname
     $        .eq.'  U').or.(resname.eq.' TS').or.(resname.eq.' GS').or.
     $        (resname.eq.' SC').or.(resname.eq.' AS')).or.
     $   (resname.eq.'  N').or.    
     $   (resname.eq.'  S').or.    
     $   (resname.eq.' +A').or.    
     $   (resname.eq.' +C').or.    
     $   (resname.eq.' +G').or.    
     $   (resname.eq.' +I').or.    
     $   (resname.eq.' +T').or.    
     $   (resname.eq.' +U').or.    
     $   (resname.eq.' PU').or.    
     $   (resname.eq.' YG').or.    
     $   (resname.eq.'1AP').or.    
     $   (resname.eq.'1MA').or.    
     $   (resname.eq.'1MG').or.    
     $   (resname.eq.'2DA').or.    
     $   (resname.eq.'2DT').or.    
     $   (resname.eq.'2MA').or.    
     $   (resname.eq.'2MG').or.    
     $   (resname.eq.'4SC').or.    
     $   (resname.eq.'4SU').or.    
     $   (resname.eq.'5IU').or.    
     $   (resname.eq.'5MC').or.    
     $   (resname.eq.'5MU').or.    
     $   (resname.eq.'5NC').or.    
     $   (resname.eq.'6MP').or.    
     $   (resname.eq.'7MG').or.    
     $   (resname.eq.'A23').or.    
     $   (resname.eq.'AD2').or.    
     $   (resname.eq.'AET').or.    
     $   (resname.eq.'AMD').or.    
     $   (resname.eq.'AMP').or.    
     $   (resname.eq.'APN').or.    
     $   (resname.eq.'ATP').or.    
     $   (resname.eq.'AZT').or.    
     $   (resname.eq.'CCC').or.    
     $   (resname.eq.'CMP').or.    
     $   (resname.eq.'CPN').or.    
     $   (resname.eq.'DAD').or.    
     $   (resname.eq.'DCT').or.    
     $   (resname.eq.'DDG').or.    
     $   (resname.eq.'DG3').or.    
     $   (resname.eq.'DHU').or.    
     $   (resname.eq.'DOC').or.    
     $   (resname.eq.'EDA').or.    
     $   (resname.eq.'G7M').or.    
     $   (resname.eq.'GDP').or.    
     $   (resname.eq.'GNP').or.    
     $   (resname.eq.'GPN').or.    
     $   (resname.eq.'GTP').or.    
     $   (resname.eq.'GUN').or.    
     $   (resname.eq.'H2U').or.    
     $   (resname.eq.'HPA').or.    
     $   (resname.eq.'IPN').or.    
     $   (resname.eq.'M2G').or.    
     $   (resname.eq.'MGT').or.    
     $   (resname.eq.'MIA').or.    
     $   (resname.eq.'OMC').or.    
     $   (resname.eq.'OMG').or.    
     $   (resname.eq.'OMU').or.    
     $   (resname.eq.'ONE').or.    
     $   (resname.eq.'P2U').or.    
     $   (resname.eq.'PGP').or.    
     $   (resname.eq.'PPU').or.    
     $   (resname.eq.'PRN').or.    
     $   (resname.eq.'PST').or.    
     $   (resname.eq.'QSI').or.    
     $   (resname.eq.'QUO').or.    
     $   (resname.eq.'RIA').or.    
     $   (resname.eq.'SAH').or.    
     $   (resname.eq.'SAM').or.    
     $   (resname.eq.'T23').or.    
     $   (resname.eq.'T6A').or.    
     $   (resname.eq.'TAF').or.    
     $   (resname.eq.'TLC').or.    
     $   (resname.eq.'TPN').or.    
     $   (resname.eq.'TSP').or.    
     $   (resname.eq.'TTP').or.    
     $   (resname.eq.'UCP').or.    
     $   (resname.eq.'VAA').or.    
     $   (resname.eq.'YYG').or.    
     $   (resname.eq.'70U').or.    
     $   (resname.eq.'12A').or.    
     $   (resname.eq.'2MU').or.    
     $   (resname.eq.'127').or.    
     $   (resname.eq.'125').or.    
     $   (resname.eq.'126').or.    
     $   (resname.eq.'MEP').or.    
     $   (resname.eq.'TLN').or.    
     $   (resname.eq.'ADP').or.    
     $   (resname.eq.'TTE').or.    
     $   (resname.eq.'PYO').or.    
     $   (resname.eq.'SUR').or.    
     $   (resname.eq.'PSD').or.    
     $   (resname.eq.'S4U').or.    
     $   (resname.eq.'CP1').or.    
     $   (resname.eq.'TP1').or.    
     $   (resname.eq.'NEA').or.    
     $   (resname.eq.'GCK').or.    
     $   (resname.eq.' CH').or.    
     $   (resname.eq.'EDC').or.    
     $   (resname.eq.'DFC').or.    
     $   (resname.eq.'DFG').or.    
     $   (resname.eq.'DRT').or.    
     $   (resname.eq.'2AR').or.    
     $   (resname.eq.'8OG').or.    
     $   (resname.eq.' IG').or.    
     $   (resname.eq.' IC').or.    
     $   (resname.eq.'IGU').or.    
     $   (resname.eq.'IMC').or.    
     $   (resname.eq.'GAO').or.    
     $   (resname.eq.'UAR').or.    
     $   (resname.eq.'CAR').or.    
     $   (resname.eq.'PPZ').or.    
     $   (resname.eq.'M1G').or.    
     $   (resname.eq.'ABR').or.    
     $   (resname.eq.'ABS').or.    
     $   (resname.eq.'S6G').or.    
     $   (resname.eq.'HEU').or.    
     $   (resname.eq.'  P').or.    
     $   (resname.eq.'DNR').or.    
     $   (resname.eq.'MCY').or.    
     $   (resname.eq.'TCP').or.    
     $   (resname.eq.'LGP').or.    
     $   (resname.eq.'GSR').or.    
     $   (resname.eq.'  X').or.    
     $   (resname.eq.'  R').or.    
     $   (resname.eq.'  Y').or.    
     $   (resname.eq.'  E').or.    
     $   (resname.eq.'GSS').or.    
     $   (resname.eq.'THX').or.    
     $   (resname.eq.'6CT').or.    
     $   (resname.eq.'TEP').or.    
     $   (resname.eq.'GN7').or.    
     $   (resname.eq.'FAG').or.    
     $   (resname.eq.'PDU').or.    
     $   (resname.eq.'MA6').or.    
     $   (resname.eq.'UMP').or.    
     $   (resname.eq.' SC').or.    
     $   (resname.eq.' GS').or.    
     $   (resname.eq.' TS').or.    
     $   (resname.eq.' AS').or.    
     $   (resname.eq.'ATD').or.    
     $   (resname.eq.'T3P').or.    
     $   (resname.eq.'5AT').or.    
     $   (resname.eq.'MMT').or.    
     $   (resname.eq.'SRA').or.    
     $   (resname.eq.'6HG').or.    
     $   (resname.eq.'6HC').or.    
     $   (resname.eq.'6HT').or.    
     $   (resname.eq.'6HA').or.    
     $   (resname.eq.'55C').or.    
     $   (resname.eq.'U8U').or.    
     $   (resname.eq.'BRO').or.    
     $   (resname.eq.'BRU').or.    
     $   (resname.eq.'5IT').or.    
     $   (resname.eq.'ADI').or.    
     $   (resname.eq.'5CM').or.    
     $   (resname.eq.'IMP').or.    
     $   (resname.eq.'THM').or.    
     $   (resname.eq.'URI').or.    
     $   (resname.eq.'AMO').or.    
     $   (resname.eq.'FHU').or.    
     $   (resname.eq.'TSB').or.    
     $   (resname.eq.'CMR').or.    
     $   (resname.eq.'RMP').or.    
     $   (resname.eq.'SMP')) then
            if((string.eq.'ATOM').and.(numero_residui_ATOM.ne.nres))then
               incremento_ATOM=incremento_ATOM+1
               numero_residui_ATOM=nres
            elseif((string.eq.'HETATM').and.(numero_residui_HETATM.ne.
     $              nres)) then
               incremento_HETATM=incremento_HETATM+1
               numero_residui_HETATM=nres
            else
               continue
            endif
         else
            continue
         endif
   
       
        
         if((atname.eq.'CA ').or.(atname.eq.'P   ').or.(atname.eq.'  P '
     $        ).or.(atname.eq.'   P').or.(atname.eq.' P  ')) then
            if(occupancy.eq.1.) then
               nca=nca+1 
                
c------------------------------------------
c            print*,nca,i
c------------------------------------------

               do k=1,3
                  pos(k,nca)=tau(k)
               enddo

               ch(nca)=chain

               anam(nca)=atname
               renam(nca)=resname

               num_res(nca)=nres !Per ciascun atomo di P (nca conta quanti ne ho) associo il suo nres cosicche poi posso trovarli con facilita'
               
               
c------------------------------------------
c            print*,anam(nca),renam(nca),(pos(k,nca),k=1,3)
c------------------------------------------
             elseif(occupancy.lt.1) then 
                useless=useless+1 !aggiorno il contatore useless invece che nca

             else
                continue
c----puo servire se trovo un modo di trovare una specie di posizione media fra tutte le posizioni alternative dell'atomo frazionario sfruttando il dato del fattore
c----temperatura.
c                  do k=1,3
c                   pos_useless(k,nca)=????
c                  enddo
c--------------------------------------------------------------------------------------------------------

             endif
          endif
 10       continue
       enddo
 11    close(12)

       numero_residui=incremento_ATOM+incremento_HETATM
c---------------------Il seguente e' un check per vedere se dal file .inp e dalla lettura del file pdb ricavo le stesse informazioni
      
       if(NUM_STEP.eq.numero_residui/2-1) then
          print*,'Il file ha superato la prova'
          condizione=.false.
       else
          differenza=(numero_residui/2-1)-NUM_STEP
          print*,'Qualche problema; nel file mancano ',differenza,
     $         ' coppie di basi'
           condizione=.true.
       endif
       







c      goto 20

c      goto 23

c      goto 25

c      goto 30

c     goto 35 !angle for (P and Ca)

c      goto 40  !dihedral

c      goto 45

c      goto 55

c      goto 60   !non bonded

c      goto 65    !secondary
            
       goto 70
      
c      goto 80

c      goto 90

c      goto 100

c      goto 110

c      goto 120
      
c      goto 130




 70   continue



      print*,incremento_ATOM,incremento_HETATM

      print*,nca
    
c      do i=1,nca
c        print*,num_res(i)
c      enddo

c----------------------------------------------------
c     First Neighbours
c----------------------------------------------------
      open(11,file='NEW_PPtipoA.dat')


      
     

      
      do i=1,nd
         
         do ii=1,nca/2
            
              
           
               
               do jj=ii+1,nca/2
               
c                    print*,'*****'
c                    print*,num_res(ii)
c                    print*,num_res(jj)
                     rr=sqrt((pos(1,ii)-pos(1,jj))**2+
     $                (pos(2,ii)-pos(2,jj))**2+(pos(3,ii)-pos(3,jj))**2)
                  
                   

c------- !MI SERVE PER SAPERE QUALE FILE CONTIENE IL DATO "SBAGLIATO" PERCHE LA DISTRIBUZIONE IN XMGRACE PUO CONTENERE UNO SPIKE IN 1
c               if((rr.gt.1-rd/2.).and.(rr.le.1+rd/2.)) then
c                  print*,filename                            
c               endif
c-------------------------------------------------------------------------------------------------------------------------------
                     if((rr.gt.r(i)-rd/2.).and.(rr.le.r(i)+rd/2.)) then

                        
                        nam1=anam(ii)
                        nam2=anam(jj)



                        if((nam1(1:2).eq.'P ').and.(nam2(1:2).eq.'P ')
     $                       .or.(nam1(2:3).eq.'P ').and.(nam2(2:3)
     $                       .eq.'P ').or.(nam1(3:4).eq.'P ').and.   !HO AGGIUNTO QUALCHE CONDIZIONE PER VIA DEI CAMBIAMENTI FATTI SUL FORMAT ECC...
     $                       (nam2(3:4).eq.'P ')) then

                        
                           if((jj.eq.ii+1).and.(num_res(jj).eq.
     $                      num_res(ii)+1).and.(ch(ii).eq.ch(jj))) then !seleziona solo i primi vicini reali cioe' 
                                                                         !1) con jj=ii+1
                                                                         !2) considerando che possono mancare dei P nella catena e quindi il buco mi darebbe una distanza fra due P adiacenti 
                                                                         !sbagliata
                                                                         !3) lungo la stessa catena
                              

                              do k=1,NUM_STEP
c                print*,num_res(ii),num_res(jj),residuo1(k),residuo3(k)
                              if((num_res(ii).eq.residuo1(k)).and.
     $                             (num_res(jj).eq.residuo3(k))) then
c                                 print*,'1-3',num_res(ii),num_res(jj)
c                                 print*,k,jj,ii,i
                                 caca(i)=caca(i)+1
c                                 print*,rr
                              else
                                 continue
                              endif
c                 print*,num_res(ii),num_res(jj),residuo2(k),residuo4(k)
                             
                             
                             enddo
                           elseif(jj.eq.ii+1) then
c                              print*,num_res(jj),num_res(ii) !serve solo come check per vedere se fa la cosa giusta
                           endif
                        
                        endif
                     endif
       
                  
               enddo

             
         enddo
c         write(11,*) r(i),caca(i)


      enddo
      
      do i=1,nd

         do ii=nca/2+1,nca
            
              
           
               
            do jj=ii+1,nca
               
c     print*,'*****'
c     print*,num_res(ii)
c     print*,num_res(jj)
               rr=sqrt((pos(1,ii)-pos(1,jj))**2+
     $              (pos(2,ii)-pos(2,jj))**2+(pos(3,ii)-pos(3,jj))**2)
               
               
               
c------- !MI SERVE PER SAPERE QUALE FILE CONTIENE IL DATO "SBAGLIATO" PERCHE LA DISTRIBUZIONE IN XMGRACE PUO CONTENERE UNO SPIKE IN 1
c     if((rr.gt.1-rd/2.).and.(rr.le.1+rd/2.)) then
c     print*,filename                            
c     endif
c-------------------------------------------------------------------------------------------------------------------------------
               if((rr.gt.r(i)-rd/2.).and.(rr.le.r(i)+rd/2.)) then
                  
                  
                  nam1=anam(ii)
                  nam2=anam(jj)
                  
                  
                  
                  if((nam1(1:2).eq.'P ').and.(nam2(1:2).eq.'P ')
     $                 .or.(nam1(2:3).eq.'P ').and.(nam2(2:3)
     $                 .eq.'P ').or.(nam1(3:4).eq.'P ').and. !HO AGGIUNTO QUALCHE CONDIZIONE PER VIA DEI CAMBIAMENTI FATTI SUL FORMAT ECC...
     $                 (nam2(3:4).eq.'P ')) then
                     
                     
                     if((jj.eq.ii+1).and.(num_res(jj).eq.
     $                    num_res(ii)+1).and.(ch(ii).eq.ch(jj))) then !seleziona solo i primi vicini reali cioe' 
                                !1) con jj=ii+1
                                !2) considerando che possono mancare dei P nella catena e quindi il buco mi darebbe una distanza fra due P adiacenti 
                                !sbagliata
                                !3) lungo la stessa catena
                        do k=1,NUM_STEP 
c                  print*,num_res(ii),num_res(jj),residuo4(k),residuo2(k)
                           if((num_res(ii).eq.residuo4(k)).and.
     $                          (num_res(jj).eq.residuo2(k))) then
c                              print*,'2-4',num_res(ii),num_res(jj)
c                              print*,k,jj,ii,i
                              caca(i)=caca(i)+1
c     print*,rr
                           else
                              continue
                           endif
                        enddo
                     elseif(jj.eq.ii+1) then
c                              print*,num_res(jj),num_res(ii) !serve solo come check per vedere se fa la cosa giusta
                     endif
                     
                  endif
               endif
       
               
            enddo

             
         enddo
         write(11,*) r(i),caca(i)
         
         

      enddo
                        
      close(11)

   
     
      close(12) !CHIUDE IL FILE PDB CORRENTE

      do i=1,2000
         residuo1(i)=-1000000
         residuo2(i)=-1000000
         residuo3(i)=-1000000
         residuo4(i)=-1000000

      enddo
      
 71   goto 5
      
c---------------------------------------------------------------------------
      stop
c---------------------------------------------------------------------------

 999   close(10)
       close(15)
       stop
 1     print*,'Manca il file percentuali'
       stop
 1000  format(a6,i5,1x,a4,1x,a3,1x,a1,i4,4x,3f8.3,2f6.3)!HO MODIFICATO PERCHE' PRIMA ERA format(a6,i5,2x,a3,1x,a3,1x,a1,i4,4x,3f8.3,2f6.3) E QUINDI HO MODIFICATO ANCHE atname e  L'IF(ATNAME....)
 2000  format(a6,i5,2x,a3,1x,a3,1x,a1,i4,4x,3f8.3,2f6.3)
 3000  format(a6,i5,2x,a2,2x,a3,1x,a1,i4,4x,3f8.3,2f6.3)
 1010  format(6x,10a5)
 1011  format(f6.3,10i5)
 1020  format(f10.5,10i5,f10.5)
 1030  format(a30)
 1040  format(4f8.3,i10)
 1050  format(3f8.3,i5)
      
 110  format(a7,16x,i3)
 120  format(a4)
 200  format(25x,a5,19x,a5)
 300  format(28x,i1,23x,i1)
 400  format(28x,i1,22x,i2)
 500  format(28x,i1,21x,i3)
 550  format(27x,i2,23x,i1)
 600  format(27x,i2,22x,i2)
 700  format(27x,i2,21x,i3)
 730  format(26x,i3,23x,i1)
 760  format(26x,i3,22x,i2)
 800  format(26x,i3,21x,i3)
      end
