!  Transition_Pathway.f90 
!
!  FUNCTIONS:    From the output of last step, generating the transition pathway (A->B) 
!

!****************************************************************************
!
!  PROGRAM: Transition_Pathway
!
!  Author: Zheng Yang, proteinyang@gmail.com 
!
!  Date:  11/04/2006
!
!**************************************************************************************** 

        program Transition_Pathway

        implicit none
        
! ---- Variables : *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Part 1 : Read in [Transition_Pathway.inp]

       Character*4    PDBid_A   
       Character*4    PDBid_B   
       
       Integer        Nitr    ! Number of total iterations
       Real           RMSDitv ! RMSD interval in the transition pathway (0:chose the smallest RMSD)
       Integer        Mstart  ! The starting MODEL number
       Character*120  Chartemp
       
! _ _ _ _ _ _ _P_D_B_ _F_I_L_E_ _R_E_A_D_I_N_G_ _V_A_R_I_A_B_L_E_S_:_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !
! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
!  Variables: Parameter
       
! Variables : ------------------------------------------------------- Part 3.1 : Scan the PDB Files  [Atom Part]
       Integer    PDBerror_A         ! The errorcode druing PDB file reading 
       Integer    MinModel_A, MinChain_A, MinResidue_A, MinAtom_A  ! The initial Min value is set to an extreme big value
       Integer    MaxModel_A, MaxChain_A, MaxResidue_A, MaxAtom_A  ! The initial Max value is set to an extreme small value
       Integer    MaxHelix_A, MaxSheet_A, MaxStrand_A,  MaxTurn_A  ! 

! Variables : ------------------------------------------------------- Part 3.2 : Read in the PDB Files  [Atom Part]
       Character, ALLOCATABLE ::    ATM_A( : , : , : , : )*4      ! The Name of atom
       Character, ALLOCATABLE ::    AANm_A( : , : , : , : )*3     ! The Name of residue
       Real, ALLOCATABLE ::         X_A( : , : , : , : )          ! Coordinates: X
       Real, ALLOCATABLE ::         Y_A( : , : , : , : )          ! Coordinates: Y
       Real, ALLOCATABLE ::         Z_A( : , : , : , : )          ! Coordinates: Z
       Real, ALLOCATABLE ::         B_A( : , : , : , : )          ! B-factors
       
       Integer, ALLOCATABLE ::      vd_A( : , : , : , : )         ! Record the existence of the ATOM :
                                                                  ! Amino Acids Atom   : vd=1
                                                                  ! Nucleic Acids Atom : vd=2

! Variables : ------------------------------------------------------- Part 3.3 : Read in the PDB Files  [Head Part]

       Character    PDBRecord_A(100,5000)*80 ! Recording the content of each line : PDBRecord(RecordType, Line#) 
       Integer      Line_A(100)              ! The line number of each RecordType : Line(RecordType)
                                                                
                                                                           
! Variables : ------------------------------------------------------- Part 3.5 : Put the coordinates in a sequence
       Integer     Site_A       ! The total nodes number
       Integer,    Allocatable :: SeqIndex_A( : , : )     ! Array SeqIndex(Site#,1/2/3/4) /Model#/Chain#/Residue#/Atom#/
       Integer,    Allocatable :: Seq_A( : , : , : , : )  ! The sequential number of a certain node at (k,j,i,s)
       
! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !       
! _ _ _ _ _ _ _P_D_B_ _F_I_L_E_ _R_E_A_D_I_N_G_ _V_A_R_I_A_B_L_E_S_:_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !
! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
       
! Variables : ------------------------------------------------------- Part 4.1 : Scan the PDB Files  [Atom Part]
       Integer    PDBerror_B         ! The errorcode druing PDB file reading 
       Integer    MinModel_B, MinChain_B, MinResidue_B, MinAtom_B  ! The initial Min value is set to an extreme big value
       Integer    MaxModel_B, MaxChain_B, MaxResidue_B, MaxAtom_B  ! The initial Max value is set to an extreme small value
       Integer    MaxHelix_B, MaxSheet_B, MaxStrand_B,  MaxTurn_B  ! 

! Variables : ------------------------------------------------------- Part 4.2 : Read in the PDB Files  [Atom Part]
       Character, ALLOCATABLE ::    ATM_B( : , : , : , : )*4      ! The Name of atom
       Character, ALLOCATABLE ::    AANm_B( : , : , : , : )*3     ! The Name of residue
       Real, ALLOCATABLE ::         X_B( : , : , : , : )          ! Coordinates: X
       Real, ALLOCATABLE ::         Y_B( : , : , : , : )          ! Coordinates: Y
       Real, ALLOCATABLE ::         Z_B( : , : , : , : )          ! Coordinates: Z
       Real, ALLOCATABLE ::         B_B( : , : , : , : )          ! B-factors
       
       Integer, ALLOCATABLE ::      vd_B( : , : , : , : )         ! Record the existence of the ATOM :
                                                                  ! Amino Acids Atom   : vd=1
                                                                  ! Nucleic Acids Atom : vd=2

! Variables : ------------------------------------------------------- Part 4.3 : Read in the PDB Files  [Head Part]

       Character    PDBRecord_B(100,5000)*80 ! Recording the content of each line : PDBRecord(RecordType, Line#) 
       Integer      Line_B(100)              ! The line number of each RecordType : Line(RecordType)
                                                                
                                                                           
! Variables : ------------------------------------------------------- Part 3.5 : Put the coordinates in a sequence
       Integer     Site_B       ! The total nodes number
       Integer,    Allocatable :: SeqIndex_B( : , : )     ! Array SeqIndex(Site#,1/2/3/4) /Model#/Chain#/Residue#/Atom#/
       Integer,    Allocatable :: Seq_B( : , : , : , : )  ! The sequential number of a certain node at (k,j,i,s)
       
! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !              

! ---- Variables :  [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] Part 6.0 : Locate the common range ....        
       
       Integer              lmtA(4,2)   ! Range of Structure A
       Integer              lmtB(4,2)   ! Range of Structure B
       Integer              lmtC(4,2)   ! Range of Common 
       
       Integer              SITE
       Integer, Allocatable :: SeqIndex( : , : ) 
       Integer, Allocatable :: Seq( : , : , : , : ) 
       
       Integer              MinModel, MinChain, MinResidue, MinAtom  
       Integer              MaxModel, MaxChain, MaxResidue, MaxAtom  
       
       Integer, Allocatable :: vd(:,:,:,:)
       
! ---- Variables :  &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*&  Part 7.0 : Projecting the data to 1D array

       Real, Allocatable :: rA(:,:) ! The coordinates of the structure A
       Real, Allocatable :: rB(:,:) ! The coordinates of the structure B              
       
! ---- Variables : <*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.  Part 8.0 : Read in the deformation vector
       
       Integer     itrt
	   Character   ACE*2 ! the string will be stored here 
       
       Real, Allocatable :: dfmVA(:,:,:) ! The deformation vector for structure A 
       Real, Allocatable :: dfmVB(:,:,:) ! The deformation vector for structure B 

! ---- Variables : <RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD> Part 9.0 : Find the minimum RMSD       
       
       Real, Allocatable :: dist(:,:)
       
       Real        distMin            ! The smallest RMSD between adjacent states
       Integer     Nstp               ! The total steps in the transition pathway
       
! ---- Variables : #=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=# Part 10.0 : Generate the structures in the transition pathway
       
       Real, Allocatable :: rn(:,:,:)

! ---- Variables : T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T Part 11.0 : Output the trajectory file        
       
! ---- Variables : E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E Part 13.0 : evaluating the system energy         
       
       Real, Allocatable :: Engy(:)
       
! ---- Variables : ------------------------------------------------------------------- 
       Integer     i, j, k, s, n, m
       Integer     ct

! ===== Body of Transition_Pathway ==================================================================================================================== ! 

! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Part 1 : Read in [Transition_Pathway.inp]        
       Open (2, File='Transition_Pathway.inp')
       
       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,'(A4)') PDBid_A
       
       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,'(A4)') PDBid_B       

       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) Nitr      ! Number of total iterations  

       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) RMSDitv   ! RMSD interval in the transition pathway (0:chose the smallest RMSD)       

       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) Mstart      ! The starting MODEL number  

       Close(2) ! File='Transition_Pathway.inp'
        
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!   *   *   *   * P * D * B *   * F * I * L * E *   * R * E * A * D * I * N * G *   * M * A * I * N *   * P * R * O * G * R * A * M *   *   *   *  !       
      
! ------------------------------------------------------------------ Part 3.1 : Scan the PDB Files  [Atom Part]
      
      PDBerror_A=100
      
      Call ScanPDB(PDBid_A, MinModel_A, MaxModel_A, MinChain_A, MaxChain_A, MinResidue_A, MaxResidue_A, MinAtom_A, MaxAtom_A, MaxHelix_A, MaxSheet_A, MaxStrand_A, MaxTurn_A, PDBerror_A)
      
      
      
! ------------------------------------------------------------------ Part 3.2 : Read in the PDB Files  [Atom Part]
      
       ALLOCATE (  ATM_A( MinModel_A:MaxModel_A, MinChain_A:MaxChain_A, MinResidue_A:MaxResidue_A, MinAtom_A:MaxAtom_A ) )
       ALLOCATE ( AANm_A( MinModel_A:MaxModel_A, MinChain_A:MaxChain_A, MinResidue_A:MaxResidue_A, MinAtom_A:MaxAtom_A ) )
       ALLOCATE (    X_A( MinModel_A:MaxModel_A, MinChain_A:MaxChain_A, MinResidue_A:MaxResidue_A, MinAtom_A:MaxAtom_A ) )
       ALLOCATE (    Y_A( MinModel_A:MaxModel_A, MinChain_A:MaxChain_A, MinResidue_A:MaxResidue_A, MinAtom_A:MaxAtom_A ) )
       ALLOCATE (    Z_A( MinModel_A:MaxModel_A, MinChain_A:MaxChain_A, MinResidue_A:MaxResidue_A, MinAtom_A:MaxAtom_A ) )
       ALLOCATE (    B_A( MinModel_A:MaxModel_A, MinChain_A:MaxChain_A, MinResidue_A:MaxResidue_A, MinAtom_A:MaxAtom_A ) )
       ALLOCATE (   vd_A( MinModel_A:MaxModel_A, MinChain_A:MaxChain_A, MinResidue_A:MaxResidue_A, MinAtom_A:MaxAtom_A ) )
      
       
       Call ReadPDB(PDBid_A, MinModel_A, MaxModel_A, MinChain_A, MaxChain_A, MinResidue_A, MaxResidue_A, MinAtom_A, MaxAtom_A, PDBerror_A, ATM_A, AANm_A, X_A, Y_A, Z_A, B_A, vd_A)

! ------------------------------------------------------------------ Part 3.3 : Read in the PDB Files  [Head Part]
       
       Call HeadPDB(PDBid_A, PDBRecord_A, Line_A)
            
       
! ------------------------------------------------------------------ Part 3.5 : Put the coordinates in a sequence
       
       Call SitePDB(PDBid_A, MinModel_A, MaxModel_A, MinChain_A, MaxChain_A, MinResidue_A, MaxResidue_A, MinAtom_A, MaxAtom_A, ATM_A, AANm_A, vd_A, Site_A)
       
       Allocate ( SeqIndex_A(1:Site_A, 1:4) )
       Allocate ( Seq_A( MinModel_A:MaxModel_A, MinChain_A:MaxChain_A, MinResidue_A:MaxResidue_A, MinAtom_A:MaxAtom_A ) )
       
       Call SequPDB(PDBid_A, MinModel_A, MaxModel_A, MinChain_A, MaxChain_A, MinResidue_A, MaxResidue_A, MinAtom_A, MaxAtom_A, ATM_A, AANm_A, vd_A, Site_A, SeqIndex_A, Seq_A)
       
       
!   *   *   *   * P * D * B *   * F * I * L * E *   * R * E * A * D * I * N * G *   * M * A * I * N *   * P * R * O * G * R * A * M *   *   *   *  !       
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!   *   *   *   * P * D * B *   * F * I * L * E *   * R * E * A * D * I * N * G *   * M * A * I * N *   * P * R * O * G * R * A * M *   *   *   *  !       
      
! ------------------------------------------------------------------ Part 4.1 : Scan the PDB Files  [Atom Part]
      
      PDBerror_B=100
      
      Call ScanPDB(PDBid_B, MinModel_B, MaxModel_B, MinChain_B, MaxChain_B, MinResidue_B, MaxResidue_B, MinAtom_B, MaxAtom_B, MaxHelix_B, MaxSheet_B, MaxStrand_B, MaxTurn_B, PDBerror_B)
      
      
      
! ------------------------------------------------------------------ Part 4.2 : Read in the PDB Files  [Atom Part]
      
       ALLOCATE (  ATM_B( MinModel_B:MaxModel_B, MinChain_B:MaxChain_B, MinResidue_B:MaxResidue_B, MinAtom_B:MaxAtom_B ) )
       ALLOCATE ( AANm_B( MinModel_B:MaxModel_B, MinChain_B:MaxChain_B, MinResidue_B:MaxResidue_B, MinAtom_B:MaxAtom_B ) )
       ALLOCATE (    X_B( MinModel_B:MaxModel_B, MinChain_B:MaxChain_B, MinResidue_B:MaxResidue_B, MinAtom_B:MaxAtom_B ) )
       ALLOCATE (    Y_B( MinModel_B:MaxModel_B, MinChain_B:MaxChain_B, MinResidue_B:MaxResidue_B, MinAtom_B:MaxAtom_B ) )
       ALLOCATE (    Z_B( MinModel_B:MaxModel_B, MinChain_B:MaxChain_B, MinResidue_B:MaxResidue_B, MinAtom_B:MaxAtom_B ) )
       ALLOCATE (    B_B( MinModel_B:MaxModel_B, MinChain_B:MaxChain_B, MinResidue_B:MaxResidue_B, MinAtom_B:MaxAtom_B ) )
       ALLOCATE (   vd_B( MinModel_B:MaxModel_B, MinChain_B:MaxChain_B, MinResidue_B:MaxResidue_B, MinAtom_B:MaxAtom_B ) )
      
       
       Call ReadPDB(PDBid_B, MinModel_B, MaxModel_B, MinChain_B, MaxChain_B, MinResidue_B, MaxResidue_B, MinAtom_B, MaxAtom_B, PDBerror_B, ATM_B, AANm_B, X_B, Y_B, Z_B, B_B, vd_B)

! ------------------------------------------------------------------ Part 4.3 : Read in the PDB Files  [Head Part]
       
       Call HeadPDB(PDBid_B, PDBRecord_B, Line_B)
            
       
! ------------------------------------------------------------------ Part 4.5 : Put the coordinates in a sequence
       
       Call SitePDB(PDBid_B, MinModel_B, MaxModel_B, MinChain_B, MaxChain_B, MinResidue_B, MaxResidue_B, MinAtom_B, MaxAtom_B, ATM_B, AANm_B, vd_B, Site_B)
       
       Allocate ( SeqIndex_B(1:Site_B, 1:4) )
       Allocate ( Seq_B( MinModel_B:MaxModel_B, MinChain_B:MaxChain_B, MinResidue_B:MaxResidue_B, MinAtom_B:MaxAtom_B ) )
       
       Call SequPDB(PDBid_B, MinModel_B, MaxModel_B, MinChain_B, MaxChain_B, MinResidue_B, MaxResidue_B, MinAtom_B, MaxAtom_B, ATM_B, AANm_B, vd_B, Site_B, SeqIndex_B, Seq_B)
       
       
!   *   *   *   * P * D * B *   * F * I * L * E *   * R * E * A * D * I * N * G *   * M * A * I * N *   * P * R * O * G * R * A * M *   *   *   *  !       
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!            

! [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] [-c-] Part 6.0 : Locate the common range ....        
              
       lmtA(1,1)=MinModel_A
       lmtA(1,2)=MaxModel_A
       lmtA(2,1)=MinChain_A
       lmtA(2,2)=MaxChain_A
       lmtA(3,1)=MinResidue_A
       lmtA(3,2)=MaxResidue_A
       lmtA(4,1)=MinAtom_A
       lmtA(4,2)=MaxAtom_A

       lmtB(1,1)=MinModel_B
       lmtB(1,2)=MaxModel_B
       lmtB(2,1)=MinChain_B
       lmtB(2,2)=MaxChain_B
       lmtB(3,1)=MinResidue_B
       lmtB(3,2)=MaxResidue_B
       lmtB(4,1)=MinAtom_B
       lmtB(4,2)=MaxAtom_B
       
       Call CommonRange ( lmtA, lmtB, lmtC, SITE_A, SITE_B, SITE )

       MinModel      =  lmtC(1,1)
       MaxModel      =  lmtC(1,2)
       MinChain      =  lmtC(2,1)
       MaxChain      =  lmtC(2,2)
       MinResidue    =  lmtC(3,1)
       MaxResidue    =  lmtC(3,2)
       MinAtom       =  lmtC(4,1)
       MaxAtom       =  lmtC(4,2)

       Allocate( vd(MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom) )
       Allocate( SeqIndex(1:SITE, 1:4) )
       Allocate( Seq(MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom) )

       Call CommonResidue ( lmtA, lmtB, lmtC, vd_A, vd_B, AANm_A, AANm_B, ATM_A, ATM_B, vd, SITE, SeqIndex, Seq )

! &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*& &*&  Part 7.0 : Projecting the data to 1D array
       
       Allocate( rA(1:SITE, 1:3) )
       Allocate( rB(1:SITE, 1:3) )
       
       Do n=1, SITE
          
          k=SeqIndex(n,1)          
          j=SeqIndex(n,2)          
          i=SeqIndex(n,3)          
          s=SeqIndex(n,4)          
          
          rA(n,1) = x_A(k,j,i,s)
          rA(n,2) = y_A(k,j,i,s)
          rA(n,3) = z_A(k,j,i,s)

          rB(n,1) = x_B(k,j,i,s)
          rB(n,2) = y_B(k,j,i,s)
          rB(n,3) = z_B(k,j,i,s)                    
          
       ENddo       
       
!       dist0 = 1.0 * rmsd(SITE, rA, rB) * Sqrt(1.0*SITE)
!       
!       Write(*,*) '                                     ' 
!       Write(*,*) 'The initial RMSD is : ', rmsd(SITE, rA, rB)
!       Write(*,*) '                                     '              
        
! <*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*>.<*> Part 8.0 : Read in the deformation vector
        
        If ( Nitr>0 ) Then
            Allocate( dfmVA(1:Nitr, 1:SITE, 1:3) )
            Allocate( dfmVB(1:Nitr, 1:SITE, 1:3) )
        Else If(Nitr==0) Then
            Allocate( dfmVA(1:1, 1:SITE, 1:3) )
            Allocate( dfmVB(1:1, 1:SITE, 1:3) )
        EndIf
        
        If( Nitr>0 ) Then 
        
            Do itrt=1, Nitr
                
                Open(800+itrt, File='['//PDBid_A//']_DfmVector_Iteration'//ACE(itrt,2)//'.xyz' )
                   
                Do n=1, SITE
                    Read(800+itrt,*) i, dfmVA(itrt,n,1), dfmVA(itrt,n,2), dfmVA(itrt,n,3)
                EndDo
                   
                Close(800+itrt) ! '['//PDBid_A//']_DfmVector_Iteration'//ACE(itrt,2)//'.xyz'
                
                ! ---------------------------------------------------------------------- ! 

                Open(850+itrt, File='['//PDBid_B//']_DfmVector_Iteration'//ACE(itrt,2)//'.xyz' )
                   
                Do n=1, SITE
                    Read(850+itrt,*) i, dfmVB(itrt,n,1), dfmVB(itrt,n,2), dfmVB(itrt,n,3)
                EndDo
                   
                Close(850+itrt) ! '['//PDBid_B//']_DfmVector_Iteration'//ACE(itrt,2)//'.xyz'
                
            
            Enddo ! itrt=1, Nitr
        
        Else If(Nitr==0) Then
        
            Do n=1, SITE
               dfmVA(1,n,1)=0.0
               dfmVA(1,n,2)=0.0
               dfmVA(1,n,3)=0.0
               dfmVB(1,n,1)=0.0
               dfmVB(1,n,2)=0.0
               dfmVB(1,n,3)=0.0
            EndDo ! n
            Nitr=1
        EndIF

! <RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD>.<RMSD> Part 9.0 : Find the minimum RMSD       
        
        Allocate( dist(0:Nitr, 1:2) )
        
        Call MinRMSD ( SITE, Nitr, rA, rB, dfmVA, dfmVB, dist, distMin ) 
        
        If(RMSDitv==0.0) RMSDitv=distMin
        
        If(RMSDitv<=1.0E-4) Then
            Write(*,'(A)') 'Error!  The RMSD interval in the transition pathway is too small ! '
            Stop 
        EndIf 
        
        ct=0
        Do k=1, Nitr
           If( dist(k,1)>=RMSDitv ) Then
              ct = ct + dist(k,1)/RMSDitv
           Else
              ct = ct +1
           EndIf
           If( dist(k,2)>=RMSDitv ) Then
              ct = ct + dist(k,2)/RMSDitv
           Else
              ct = ct +1
           EndIf
        Enddo
        
        If( dist(0,2)>=RMSDitv ) then
            ct = ct + dist(0,2)/RMSDitv +1 
        Else
            ct = ct +2
        EndIf
        
        Nstp=ct
        
        Write(*,'(A,I6)') 'The steps in the transition pathway will be : ', Nstp
        Write(*,'(A)')'                             '
        
        Write(*,'(A)') 'Do you want to continue? (Y/N)' 
        Read(*,'(A1)') Chartemp(1:1)
        
        If( Chartemp(1:1)=='N' .OR. Chartemp(1:1)=='n' ) Goto 777

! #=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=# Part 10.0 : Generate the structures in the transition pathway      
        
        Allocate ( rn(1:Nstp, 1:SITE, 1:3) ) ! The coordinates of the strucutre in the transition pathway
        
        Call InterStru ( SITE, Nitr, Nstp, RMSDitv, rA, rB, dfmVA, dfmVB, dist, rn ) 

! T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T-T Part 11.0 : Output the trajectory file        
        
        Call TrajectoryWriting ( PDBid_A, PDBid_B, lmtA, SITE, Nstp, PDBRecord_A, Line_A, ATM_A, AANm_A, X_A, Y_A, Z_A, vd_A, Seq_A, SeqIndex_A, rn, Mstart ) 
        
! P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P*P Part 12.0 : Generate a PyMOL script

        Call TransPyMOL ( PDBid_A, PDBid_B, Nstp ) 
        
! E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E+E Part 13.0 : evaluating the system energy         
        
        Allocate ( Engy(1:Nstp) )
        
        Call EnergyEvaluating ( PDBid_A, PDBid_B, SITE, Nstp, Nitr, rn, Engy )
        
        Call TopEngStruc ( PDBid_A, PDBid_B, lmtA, SITE, Nstp, PDBRecord_A, Line_A, ATM_A, AANm_A, X_A, Y_A, Z_A, vd_A, Seq_A, SeqIndex_A, rn, Engy )
        
! ===== E = N = D ===================================================================================================================================== !
777     Stop        
        end program Transition_Pathway



























