!  Intermediate_States_Generator.f90 
!
!  FUNCTIONS:
!  Intermediate_States_Generator      - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Intermediate_States_Generator
!
!  Author: Zheng Yang, proteinyang@gmail.com 
!
!  Date:  8/26/2006
!
!**************************************************************************************** 

       program Intermediate_States_Generator

       implicit none


! ---- Variables : *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Part 1 : Read in [Intermediate_States_Generator.inp]

       Character*4    PDBid_A   
       Character*4    PDBid_B   
       
       Integer        MEIG    ! The total eigen modes solved by ANM
       
       Real           Dmntmd  ! The threshold of the domint modes
       Integer        itrt    ! The iteration number
       Integer        Nstep   ! The trajectory steps 
       Real           FmaxA   ! The cumulated correlation threshold
       Real           FmaxB   ! The cumulated correlation threshold
       Real           scalf   ! The scaling factor s in Eq.4
       Real           rmsdd0  ! The initial RMSD between ending strcutures
       
       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 : *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- Part 5.0 : Read in the Eigenvalues & Eigenvectors of Structure A

       Real, ALLOCATABLE :: vx_A(:,:), vy_A(:,:), vz_A(:,:)   !  Eigenvectors 
       Real, ALLOCATABLE :: vx_B(:,:), vy_B(:,:), vz_B(:,:)   !  Eigenvectors 

! ---- 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
              
       Real, Allocatable :: vA(:,:,:) ! The eigenvectors of the structure A              
       Real, Allocatable :: vB(:,:,:) ! The eigenvectors of the structure B  
       
       Integer              na, nb            
       
       Real                 rmsd
       Real                 dist0  ! the magnitude of the initial displacement vector

! ---- Variables : ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ Part 7.9 : Generate the step sequence 
       
       Real,          Allocatable :: Step_A( : , : ) ! The step sequence, amplitude of the eigenvector
       Real,          Allocatable :: Step_B( : , : ) ! The step sequence, amplitude of the eigenvector

       Integer        Ncm_A            ! The total combined modes for structure A
       Integer        Ncm_B            ! The total combined modes for structure B
       
       Integer ::     ErrCode = 100    ! Error code, success = 100
       
! ---- Variables :  +=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+| Part 8.0 : Generating the intermediate states
       Character            ACE*2            ! Function: Return a string for a given integer
       Real, Allocatable :: rmsd_AB(:,:)
       Integer              mtA, mtB                  ! The minimum RMSD appears at the 'mtX' step
       Real                 rmsd_fnal        ! The rmsd between ending structures ( after scaling ) 
  
! ---- Variables :  +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+  Part 9.0 : Calculate the angle between the combined eigenvector and the displacement vector


! ---- Variables : ------------------------------------------------------------------- 
       Integer              i, j, k, s, n, m

! ==== Body of Intermediate_States_Generator ========================================================================================================= ! 


! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Part 1 : Read in [Intermediate_States_Generator.inp]

       Open (2, File='Intermediate_States_Generator.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,*) MEIG         !  The total eigen modes solved by ANM 
              
       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) Dmntmd       !  The threshold of the domint modes 

       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) itrt         !  The iteration number       
       
       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) Nstep        ! The trajectory steps 
       
       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) FmaxA       ! The cumulated correlation threshold  
       Read(2,*) FmaxB       ! The cumulated correlation threshold  
       
       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) scalf        ! The scaling factor s

       Read(2,*) CharTemp
       Read(2,*) CharTemp
       Read(2,*) rmsdd0        ! The initial RMSD between ending strcutures
       
       
       Close(2) ! Intermediate_States_Generator.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, itrt)
      
      
      
! ------------------------------------------------------------------ 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, itrt)

! ------------------------------------------------------------------ Part 3.3 : Read in the PDB Files  [Head Part]
       
       Call HeadPDB(PDBid_A, PDBRecord_A, Line_A, itrt)
            
       
! ------------------------------------------------------------------ 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, itrt)
      
      
      
! ------------------------------------------------------------------ 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, itrt)

! ------------------------------------------------------------------ Part 4.3 : Read in the PDB Files  [Head Part]
       
       Call HeadPDB(PDBid_B, PDBRecord_B, Line_B, itrt)
            
       
! ------------------------------------------------------------------ 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 *   *   *   *  !       
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!            
       
       
! *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- Part 5.1 : Read in the Eigenvalues & Eigenvectors of Structure A
       Allocate( vx_A(1:SITE_A, 1:MEIG), vy_A(1:SITE_A, 1:MEIG), vz_A(1:SITE_A, 1:MEIG) )
       
       Call Readeigen (PDBid_A, SITE_A, MEIG, vx_A, vy_A, vz_A)

! *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- Part 5.2 : Read in the Eigenvalues & Eigenvectors of Structure B
       Allocate( vx_B(1:SITE_B, 1:MEIG), vy_B(1:SITE_B, 1:MEIG), vz_B(1:SITE_B, 1:MEIG) )
       
       Call Readeigen (PDBid_B, SITE_B, MEIG, vx_B, vy_B, vz_B)
       

! [-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) )
       
       Allocate( vA(1:SITE, 1:MEIG, 1:3) )
       Allocate( vB(1:SITE, 1:MEIG, 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)          
          
          na=Seq_A(k,j,i,s)
          nb=Seq_B(k,j,i,s)
          
          Do m=1,MEIG
             vA(n,m,1) = vx_A(na,m)
             vA(n,m,2) = vy_A(na,m)
             vA(n,m,3) = vz_A(na,m)
             
             vB(n,m,1) = vx_B(nb,m)
             vB(n,m,2) = vy_B(nb,m)
             vB(n,m,3) = vz_B(nb,m)             
          Enddo
          
       ENddo       
       
       dist0 = 1.0 * rmsd(SITE, rA, rB) * Sqrt(1.0*SITE)
       
       Write(*,*) '                                     ' 
       Write(*,*) 'The initial RMSD is : ', rmsd(SITE, rA, rB)
       Write(*,*) '                                     ' 

! +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* +|* Part 7.5 : Calculate the Cosine Correlation

       Call CosCorrelation ( PDBid_A, PDBid_B, SITE, MEIG, itrt, rA, rB, vA, Dmntmd  )
       Call CosCorrelation ( PDBid_B, PDBid_A, SITE, MEIG, itrt, rB, rA, vB, Dmntmd  )

! ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ ^=^ Part 7.9 : Generate the step sequence 
       
       Allocate( Step_A(1:Nstep, 1:MEIG)  )
       Allocate( Step_B(1:Nstep, 1:MEIG)  )
       
!       Call StepSeq( PDBid_A, PDBid_B, MEIG, Nstep, itrt, Step_A, Step_B)       
       
       Call StepSeq_CumCos( PDBid_A, PDBid_B, MEIG, Nstep, itrt, FmaxA, FmaxB, Step_A, Step_B, Dmntmd, Ncm_A, Ncm_B, ErrCode )      
       
!       Call CumCos_CosDV (  PDBid_A, PDBid_B, SITE, MEIG, itrt, dist0, rA, rB, vA, vB ) 
       
! +=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+|+=+| Part 8.0 : Generating the intermediate states
       If( ErrCode == 100 ) Then
        
           Allocate( rmsd_AB(1:Nstep, 1:Nstep) ) 

           
           Call 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 ) 

           Call WriteMoviePDB ( MEIG, Nstep, dist0, itrt,   PDBid_A, SITE_A, lmtA, Step_A, PDBRecord_A, Line_A, ATM_A, AANm_A, X_A, Y_A, Z_A, vd_A, vx_A, vy_A, vz_A, Seq_A, SeqIndex_A, mtA )
           Call WriteMoviePDB ( MEIG, Nstep, dist0, itrt,   PDBid_B, SITE_B, lmtB, Step_B, PDBRecord_B, Line_B, ATM_B, AANm_B, X_B, Y_B, Z_B, vd_B, vx_B, vy_B, vz_B, Seq_B, SeqIndex_B, mtB )
           
       EndIf ! ErrCode = 100
       
! +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+ +_+  Part 9.0 : Calculate the angle between the combined eigenvector and the displacement vector
       
!       Call AngleDV ( PDBid_A, PDBid_B, SITE, MEIG, Nstep, itrt, rA, rB, vA, Step_A )
!       Call AngleDV ( PDBid_B, PDBid_A, SITE, MEIG, Nstep, itrt, rB, rA, vB, Step_B )
       
       
       Call inpRewrite( PDBid_A, PDBid_B, Ncm_A, Ncm_B, MEIG, Dmntmd, itrt, Nstep, FmaxA, FmaxB, scalf, SITE, ErrCode, rmsdd0, rmsd_fnal )
       
! ==== E == N == D ================================================================================================================================================================================================ !

       end program Intermediate_States_Generator

