!      %***********************************************************************%
!      *             Subroutine of Main, generate the step sequence            *
!      %***********************************************************************%

       Subroutine StepSeq_CumCos ( PDBid_A, PDBid_B, MEIG, Nstep, itrt, FmaxA, FmaxB, Step_A, Step_B, Dmntmd, Ncm_A, Ncm_B, ErrCode  ) 
       
       Implicit None
       
       
! ---- Variables : ---------------------------------------------------------------------------------- !       

       Character*4    PDBid_A   
       Character*4    PDBid_B        
       
       Integer        MEIG
       Integer        itrt    ! The iteration number
       Integer        Nstep
       Real           FmaxA   ! The cumulated correlation threshold
       Real           FmaxB   ! The cumulated correlation threshold
       

       Real           Step_A(1:Nstep, 1:MEIG) ! The step sequence, amplitude of the eigenvector
       Real           Step_B(1:Nstep, 1:MEIG) ! The step sequence, amplitude of the eigenvector
       
       Real           Dmntmd
       
       Integer        Ncm_A            ! The total combined modes for structure A
       Integer        Ncm_B            ! The total combined modes for structure B
       
       Integer        ErrCode          ! Error code
       
! ---- Variables :  #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# Read in [Correlation.dat]       
       
       Character*120  Chartemp
       
       Real           CorCos_A(MEIG), ACorCos_A(MEIG), CumCos_A(MEIG)
       Real           CorCos_B(MEIG), ACorCos_B(MEIG), CumCos_B(MEIG)
       
! ---- Variables :  @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ Generating the step sequence
       
       Integer        status
       
       
       Real           StepDiv

       Character      ACE*2            ! Function: Return a string for a given integer
       
! ---- Variables : ---------------------------------------------------------------------------------- !       
       Integer        n, i, j, k, ct
       
! ==== Body of Subroutine StepSeq_CumCos ============================================================ !        
       
       
! #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# #=# Read in [Correlation.dat]       

       Open (10, File=PDBid_A//' to '//PDBid_B//'_Cos  Correlation_Iterate'//ACE(itrt,2)//'.dat')
       Open (20, File=PDBid_B//' to '//PDBid_A//'_Cos  Correlation_Iterate'//ACE(itrt,2)//'.dat')
       
       Read(10,*) CharTemp
       Read(20,*) CharTemp
       
       Do n=1, MEIG
          Read(10,*) i, CorCos_A(n), Acorcos_A(n), CumCos_A(n)
          Read(20,*) i, CorCos_B(n), Acorcos_B(n), CumCos_B(n)
       Enddo
       
       
       Close(10) ! PDBid_A//' to '//PDBid_B//'_Cos  Correlation_Iterate'//ACE(itrt,2)//'.dat'
       Close(20) ! PDBid_B//' to '//PDBid_A//'_Cos  Correlation_Iterate'//ACE(itrt,2)//'.dat'
       
! @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ @^@ Generating the step sequence
       Ncm_A=0
       j=0
       status=1
       Do While (status==1)
          j=j+1
          If( CumCos_A(j)<FmaxA ) Then
              Ncm_A = Ncm_A + 1
          Else 
              Ncm_A=j
              status=0
          EndIf !  ( CumCos_A(j)<FmaxA ) 
          If(j>=MEIG) status=0   
       Enddo  ! While (status==1)
       
       ! - - - - - - - - - - - - - - - - - - !
       
       Ncm_B=0
       j=0
       status=1
       Do While (status==1)
          j=j+1
          If( CumCos_B(j)<FmaxB ) Then
              Ncm_B = Ncm_B + 1
          Else 
              Ncm_B=j
              status=0
          EndIf !  ( CumCos_B(j)<FmaxB ) 
          If(j>=MEIG) status=0   
       Enddo  ! While (status==1)
       
       ! - - - - - - - - - - - - - - - - - - !
       
       Write(*,'(A)') '            '
       Write(*,'(A)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
       
       If( Ncm_A >= MEIG )  Then
             Write(*,'(A)') 'Warning!!! You need more modes to get the threshold! [A]'
             ErrCode = 101
       EndIf
       
       If( Ncm_B >= MEIG )  Then
             Write(*,'(A)') 'Warning!!! You need more modes to get the threshold! [B]'
             ErrCode = 101
       EndIf
       
       Write(*,'(A,I)') ' For structure A, the combined modes number is : ', Ncm_A
       Write(*,'(A,I)') ' For structure B, the combined modes number is : ', Ncm_B
       
       Write(*,'(A)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
       Write(*,'(A)') '            '

! +++++++++++++++++++++++++++++++++++++++++ Re-write the Cos Correlation file

       Open (10, File=PDBid_A//' to '//PDBid_B//'_Cos  Correlation new_Iterate'//ACE(itrt,2)//'.dat')
       Open (20, File=PDBid_B//' to '//PDBid_A//'_Cos  Correlation new_Iterate'//ACE(itrt,2)//'.dat')
       
       Write(10,'(A)') '       Index    Correlation    Cosine_Correlation  Cumulated_Correlation '
       Write(20,'(A)') '       Index    Correlation    Cosine_Correlation  Cumulated_Correlation '
       
       Do n=1, MEIG
          if ( n <= Ncm_A ) Then
              Write(10,*) n, CorCos_A(n), 0.0, Acorcos_A(n), CumCos_A(n)
          Else
              Write(10,*) n, CorCos_A(n), Acorcos_A(n), 0.0, CumCos_A(n)
          EndIf
       Enddo
       
       Do n=1, MEIG
          if ( n <= Ncm_B ) Then
              Write(20,*) n, CorCos_B(n), 0.0, Acorcos_B(n), CumCos_B(n)
          else
              Write(20,*) n, CorCos_B(n), Acorcos_B(n), 0.0, CumCos_B(n)
          Endif
       Enddo
       
       Close(10) ! PDBid_A//' to '//PDBid_B//'_Cos  Correlation_Iterate'//ACE(itrt,2)//'.dat'
       Close(20) ! PDBid_B//' to '//PDBid_A//'_Cos  Correlation_Iterate'//ACE(itrt,2)//'.dat'

! [=]+[=]+[=]+[=]+[=]+[=]+[=]+[=]+[=]+[=]+[=]+[=]+[=] _Dominant modes output
       
!       Open (401, File=PDBid_A//' to '//PDBid_B//'_Dominant Modes_Iterate'//ACE(itrt,2)//'.dat')
!       Open (402, File=PDBid_B//' to '//PDBid_A//'_Dominant Modes_Iterate'//ACE(itrt,2)//'.dat')
       
!       Write(401,'(A)') ' Dominant_mode    |Cosine_Correlation|    Cosine_Correlation '
!       Write(402,'(A)') ' Dominant_mode    |Cosine_Correlation|    Cosine_Correlation '
       
       ct = 0
       Do n=1, Ncm_A
          If ( Acorcos_A(n) >= Dmntmd ) Then 
              ct = ct + 1
!              Write(401,*) n, Acorcos_A(n), CorCos_A(n)       
!+++++++              
!          Else
!              CorCos_A(n) = 0.0   
!+++++++              
          EndIf
       Enddo
!       Write(401,'(A)') '----------------------------------------------------------------------'
!       Write(401,'(A,I)') 'Total dominant modes: ', ct
       

       ct = 0
       Do n=1, Ncm_B
          If ( Acorcos_B(n) >= Dmntmd ) Then 
              ct = ct + 1
!              Write(402,*) n, Acorcos_B(n), CorCos_B(n)       
!+++++++              
!          Else
!              CorCos_B(n) = 0.0   
!+++++++                                   
          EndIf
       Enddo
!       Write(402,'(A)') '----------------------------------------------------------------------'
!       Write(402,'(A,I)') 'Total dominant modes: ', ct
       
!       Close(401)
!       Close(402)
                            
       ! = * = * = * = * = * = * = * = * = * = * = * = * = * = * = * = * = * = * = * = * = !
       
       
       Do i=1, Nstep
           Do j=1, MEIG
               If( j<=Ncm_A ) Then 
                   Step_A(i,j) = ( StepDiv(i, Nstep, 0, 1.0) ) * CorCos_A(j)
               Else 
                   Step_A(i,j) = 0.0 
               EndIf       
           Enddo ! j
       Enddo ! i
       
       ! - - - - - - - - - - - - - - - - - - !
       
       Do i=1, Nstep
           Do j=1, MEIG
               If( j<=Ncm_B ) Then 
                   Step_B(i,j) = ( StepDiv(i, Nstep, 0, 1.0) ) * CorCos_B(j)
               Else 
                   Step_B(i,j) = 0.0 
               EndIf       
           Enddo ! j
       Enddo ! i
       
       ! + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + ! 
       
!       Open(31, File=PDBid_A//' Step Sequence Check_Iterate'//ACE(itrt,2)//'.dat')
!       Open(32, File=PDBid_B//' Step Sequence Check_Iterate'//ACE(itrt,2)//'.dat')
!
!       Write(31,'(A)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
!       Write(31,'(A,I)') ' For structure A, the combined modes number is : ', Ncm_A
!       Write(31,'(A,I)') ' For structure B, the combined modes number is : ', Ncm_B
!       Write(31,'(A)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
!       Write(31,'(A)') '            '
!
!       Write(32,'(A)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
!       Write(32,'(A,I)') ' For structure A, the combined modes number is : ', Ncm_A
!       Write(32,'(A,I)') ' For structure B, the combined modes number is : ', Ncm_B
!       Write(32,'(A)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
!       Write(32,'(A)') '            '
!       
!       Do j=1, MEIG
!         
!         Write(31,'(I6,3x,11(F7.4,3x))') j, (Step_A(i,j), i=1, Nstep)
!         Write(32,'(I6,3x,11(F7.4,3x))') j, (Step_B(i,j), i=1, Nstep)
!       
!       Enddo
!
!       
!              
!       Close(31)
!       CLose(32)
       
! =====E==N===D====================================================================================== !
       Return
       End           
