!      %***************************************************************************************%
!      *          Subroutine of Main, calculating the RMSD between intermediate states         *
!      %***************************************************************************************%

       Subroutine  InterRMSD ( PDBid_A, PDBid_B, SITE, MEIG, Nstep, itrt, dist0, rA, rB, vA, vB, Step_A, Step_B, rmsd_AB, mtA, mtB, scalf, Ncm_A, Ncm_B, rmsd_fnal, FmaxA, FmaxB )

       Implicit None
       
! ---- Variables : ---------------------------------------------------------------------------------- !       
       CHARACTER   PDBid_A*4
       CHARACTER   PDBid_B*4
       
       Integer     SITE
       Integer     MEIG
       Integer     Nstep
       
       Integer     itrt    ! The iteration number
       
       Real        dist0

       Real        rA(1:SITE, 1:3) 
       Real        rB(1:SITE, 1:3) 
       
       Real        vA(1:SITE, 1:MEIG, 1:3) 
       Real        vB(1:SITE, 1:MEIG, 1:3) 
       
       Real        Step_A(1:Nstep, 1:MEIG)
       Real        Step_B(1:Nstep, 1:MEIG)
       
       Real        rmsd_AB(1:Nstep, 1:Nstep)
       
       Integer     mtA, mtB                  ! The minimum RMSD appears at the 'mtX' step
       
       Real        scalf

       Integer        Ncm_A            ! The total combined modes for structure A
       Integer        Ncm_B            ! The total combined modes for structure B

       Real           FmaxA   ! The cumulated correlation threshold
       Real           FmaxB   ! The cumulated correlation threshold
       
! ---- Variables : ---------------------------------------------------------------------------------- !       
       
       Character   ACE*2            ! Function: Return a string for a given integer
       
       Real        min
       Real        StepDiv
       
       Integer     ta, tb
       Real        sfa, sfb, sf
       
       Real        rtA(1:SITE, 1:3)
       Real        rtB(1:SITE, 1:3)
       
       Real        rmsdA_t_o(1:Nstep) ! RMSD_to_Original
       Real        rmsdA_t_t(1:Nstep) ! RMSD_to_Target

       Real        rmsdB_t_o(1:Nstep) ! RMSD_to_Original
       Real        rmsdB_t_t(1:Nstep) ! RMSD_to_Target
       
       Real        RMSD, rmsd_fnal
       
       logical     logEXIT ! the existence of the file [PDBA-PDBB]_sum.log
       
       
       Integer     n, m
       
! ==== Body of Subroutine  InterRMSD ================================================================ ! 

!       Open(80, File='['//PDBid_A//'-'//PDBid_B//'] Inter-States RMSD_Iterate'//ACE(itrt,2)//'.dat' )       

       min=dist0/Sqrt(1.0*SITE)
       mtA=1
       mtB=1
       
       Do ta=1, Nstep
       Do tb=1, Nstep
            
            Do n=1, SITE
            
                rtA(n,1) = rA(n,1)
                rtA(n,2) = rA(n,2)
                rtA(n,3) = rA(n,3)
                Do m=1, MEIG
                   rtA(n,1)  =  rtA(n,1)   +   dist0 * Step_A(ta,m) * vA(n,m,1)
                   rtA(n,2)  =  rtA(n,2)   +   dist0 * Step_A(ta,m) * vA(n,m,2)
                   rtA(n,3)  =  rtA(n,3)   +   dist0 * Step_A(ta,m) * vA(n,m,3)
                Enddo
                
                rtB(n,1) = rB(n,1)
                rtB(n,2) = rB(n,2)
                rtB(n,3) = rB(n,3)
                Do m=1, MEIG
                   rtB(n,1)  =  rtB(n,1)   +   dist0 * Step_B(tb,m) * vB(n,m,1)
                   rtB(n,2)  =  rtB(n,2)   +   dist0 * Step_B(tb,m) * vB(n,m,2)
                   rtB(n,3)  =  rtB(n,3)   +   dist0 * Step_B(tb,m) * vB(n,m,3)
                Enddo
                
            Enddo
            
            rmsd_AB(ta,tb) = rmsd( SITE, rtA, rtB )
            
            sfa = StepDiv(ta, Nstep, 0, 1.0)
            sfb = StepDiv(tb, Nstep, 0, 1.0)
            
!            Write(* ,'(2(I3,3x),2(F6.3,3x),3x,F8.2)') ta, tb, sfa, sfb, rmsd_AB(ta,tb)
!            Write(80,'(2(I3,3x),2(F6.3,3x),3x,F8.2)') ta, tb, sfa, sfb, rmsd_AB(ta,tb)
            
            If(rmsd_AB(ta,tb)<min) Then
                min=rmsd_AB(ta,tb)
                mtA=ta
                mtB=tb
            EndIf
       
       Enddo
       Enddo

       Write(*,'(A)') '----------------------------------------------------------------'
       Write(*,'(A,I3,3x,I3,3x,A,F5.2,A)') 'The smallest rmsd appears at : ', mtA, mtB, ' [RMSD=', rmsd_AB(mtA, mtB),']'
       Write(*,'(A)') '         '

! ======================== !        
       mtA = mtA * scalf
       mtB = mtB * scalf
       
       If(mtA==0) mtA=1
       If(mtB==0) mtB=1
       
       rmsd_fnal = rmsd_AB(mtA, mtB) 
! ======================== !

       Write(*,*) 'But we are going to use : ', mtA, mtB, ' [RMSD=', rmsd_fnal,']'
       Write(*,'(A)') '         '
       
       
!       Close(80) ! File='['//PDBid_A//'-'//PDBid_B//'] Inter-States RMSD_Iterate'//ACE(itrt,2)//'.dat' 
! _=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=_=      
 
       Open(81, File=PDBid_A//'_RMSD Evolving_Iterate'//ACE(itrt,2)//'.dat' )
       Open(82, File=PDBid_B//'_RMSD Evolving_Iterate'//ACE(itrt,2)//'.dat' )
       
       Write(81,'(A)') '  Step_Index   Scaling  RMSD_to_Original   RMSD_to_Target  '
       Write(82,'(A)') '  Step_Index   Scaling  RMSD_to_Original   RMSD_to_Target  '  
       
       Do ta=1, Nstep
       
          Do n=1,SITE
            rtA(n,1) = rA(n,1)
            rtA(n,2) = rA(n,2)
            rtA(n,3) = rA(n,3)
            Do m=1, MEIG
                rtA(n,1)  =  rtA(n,1)   +   dist0 * Step_A(ta,m) * vA(n,m,1)
                rtA(n,2)  =  rtA(n,2)   +   dist0 * Step_A(ta,m) * vA(n,m,2)
                rtA(n,3)  =  rtA(n,3)   +   dist0 * Step_A(ta,m) * vA(n,m,3)
            Enddo
          Enddo
          
          sfa=StepDiv(ta, Nstep, 0, 1.0)
          
          rmsdA_t_o(ta) = rmsd(SITE, rtA, rA)
          rmsdA_t_t(ta) = rmsd(SITE, rtA, rB) 
          
          Write(81,'(I8,6x,F6.3,6x,F8.3,11x,F8.3)') ta, sfa, rmsdA_t_o(ta), rmsdA_t_t(ta)          
           
       Enddo  ! ta
       
       Write(81,'(A)') '-----------------------------------------------'
       Write(81,'(A,I3,A,F5.2,A)') 'The final RMSD difference after the step :', mtA, '    [RMSD=', rmsd_fnal, ']'
       Write(81,'(A)') '-----------------------------------------------'
       Write(81,'(A)') '                                               '
       
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -    
    
       Do tb=1, Nstep
       
          Do n=1,SITE
            rtB(n,1) = rB(n,1)
            rtB(n,2) = rB(n,2)
            rtB(n,3) = rB(n,3)
            Do m=1, MEIG
                rtB(n,1)  =  rtB(n,1)   +   dist0 * Step_B(tb,m) * vB(n,m,1)
                rtB(n,2)  =  rtB(n,2)   +   dist0 * Step_B(tb,m) * vB(n,m,2)
                rtB(n,3)  =  rtB(n,3)   +   dist0 * Step_B(tb,m) * vB(n,m,3)
            Enddo
          Enddo
          
          sfb=StepDiv(tb, Nstep, 0, 1.0)

          rmsdB_t_o(tb) = rmsd(SITE, rtB, rB)
          rmsdB_t_t(tb) = rmsd(SITE, rtB, rA)           
          
          Write(82,'(I8,6x,F6.3,6x,F8.3,11x,F8.3)') tb, sfb, rmsdB_t_o(tb), rmsdB_t_t(tb)
                    
       Enddo  ! tb
       
       Write(82,'(A)') '-----------------------------------------------'
       Write(82,'(A,I3,A,F5.2,A)') 'The final RMSD difference after the step :', mtB, '    [RMSD=', rmsd_fnal, ']'
       Write(82,'(A)') '-----------------------------------------------'

              
       Close(81) ! File=PDBid_A//' RMSD Evolving.dat' 
       Close(82) ! File=PDBid_B//' RMSD Evolving.dat' 
       
       Write(*,*) ' '
       Write(*,*) 'The closest distance to the target:', PDBid_B , rmsdA_t_t(Nstep)
       Write(*,*) 'The closest distance to the target:', PDBid_A , rmsdB_t_t(Nstep)
       Write(*,*) 'The smallest rmsd if deform simultaneously:  ', rmsd_fnal
       Write(*,*) ' '

! ------------------------------------------- Write the log file
       INQUIRE ( File='['//PDBid_A//'-'//PDBid_B//']_sum.log', EXIST=logEXIT )

       Open(880, File= '['//PDBid_A//'-'//PDBid_B//']_sum.log', ACCESS = 'APPEND' )
       
       If( logEXIT == .False. ) Then
            Write(880,'(A)') '   k   m_A(k)    m_B(k)     RMSD(k)     scale_A   scale_B   |v_A(k)|  |v_B(k)|   Fmin_A    Fmin_B'
            write(880, '(I4,2(I8,2x),F10.3)') 0, 0, 0, rmsd_AB(1,1)
       endif

       sfa = StepDiv(mtA, Nstep, 0, 1.0)
       sfb = StepDiv(mtB, Nstep, 0, 1.0)
       
!       If(mtA==1) Ncm_A = 0
!       If(mtB==1) Ncm_B = 0
       
       Write(880,'(I4,2(I8,2x),F10.3,5x,6(F8.2,2x))') itrt, Ncm_A, Ncm_B, rmsd_fnal, sfa, sfb, rmsdA_t_o(mtA), rmsdB_t_o(mtB), FmaxA, FmaxB  
       
       Close(880)      

! =====E==N===D====================================================================================== !
       Return
       End                          