        !************************************************************************
        !         A Subroutine of Main, Generating the Hessian Matrix           *
        !************************************************************************
       
        Subroutine Hessian ( PDBid, rcut, gama )
       
        Implicit None
       
! ----- Variables : -----------------------------------------------------------------------!
        Character   PDBid*4
        Real        rcut      ! cutoff distance (in A)
        Real        gama   ! the force constant, gama

! _ _ _ _ _ _ _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 0.1 : Scan the PDB Files  [Atom Part]
       Integer    PDBerror         ! The errorcode druing PDB file reading 
       Integer    MinModel, MinChain, MinResidue, MinAtom  ! The initial Min value is set to an extreme big value
       Integer    MaxModel, MaxChain, MaxResidue, MaxAtom  ! The initial Max value is set to an extreme small value
       Integer    MaxHelix, MaxSheet, MaxStrand,  MaxTurn  ! 

! Variables : ------------------------------------------------------- Part 0.2 : Read in the PDB Files  [Atom Part]
       Character, ALLOCATABLE ::    ATM( : , : , : , : )*4      ! The Name of atom
       Character, ALLOCATABLE ::    AANm( : , : , : , : )*3     ! The Name of residue
       Real, ALLOCATABLE ::         X( : , : , : , : )          ! Coordinates: X
       Real, ALLOCATABLE ::         Y( : , : , : , : )          ! Coordinates: Y
       Real, ALLOCATABLE ::         Z( : , : , : , : )          ! Coordinates: Z
       Real, ALLOCATABLE ::         B( : , : , : , : )          ! B-factors
       
       Integer, ALLOCATABLE ::      vd( : , : , : , : )         ! Record the existence of the ATOM :
                                                                ! Amino Acids Atom   : vd=1
                                                                ! Nucleic Acids Atom : vd=2
                                                                
! Variables : ------------------------------------------------------- Part 0.3 : Read in the PDB Files  [Head Part]

       Character    PDBRecord(100,5000)*80 ! Recording the content of each line : PDBRecord(RecordType, Line#) 
       Integer      Line(100)             ! The line number of each RecordType : Line(RecordType)
           
! Variables : ------------------------------------------------------- Part 0.5 : Put the coordinates in a sequence
       Integer     Site       ! The total nodes number
       Integer,    Allocatable :: SeqIndex( : , : )     ! Array SeqIndex(Site#,1/2/3/4) /Model#/Chain#/Residue#/Atom#/
       Integer,    Allocatable :: Seq( : , : , : , : )  ! The sequential number of a certain node at (k,j,i,s)
       
! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !       
       
       Real(8),   Allocatable :: xx(:), yy(:), zz(:)    ! Coordinates in 1-D
       
! Part 4 --------------------------------------------------------------------------- Calculate the Hessian

       Real(8)     dn(3,3), hess
       Real(8)     bx, by, bz
       Real(8)     rij
       Integer(8)  nonzero
       Real(8)     spareness
       Real(8)     avgbond
       
       Real(8)     rijsum

       Logical     INPexist 
! ---- Variables : ---------------------------------------------------------------------------------- !       
       Integer     i, j, k, s, l, n, ctr, t, ct
        
! ===== Body of Subroutine Hessian ========================================================!

       Open(9,File=PDBid//'_Hessian.log')


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!   *   *   *   * 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 0.1 : Scan the PDB Files  [Atom Part]
      
      PDBerror=100
      
      Call ScanPDB(PDBid, MinModel, MaxModel, MinChain, MaxChain, MinResidue, MaxResidue, MinAtom, MaxAtom, MaxHelix, MaxSheet, MaxStrand, MaxTurn, PDBerror)
      
      
      
! ------------------------------------------------------------------ Part 0.2 : Read in the PDB Files  [Atom Part]
      
       ALLOCATE (  ATM( MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom ) )
       ALLOCATE ( AANm( MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom ) )
       ALLOCATE (    X( MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom ) )
       ALLOCATE (    Y( MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom ) )
       ALLOCATE (    Z( MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom ) )
       ALLOCATE (    B( MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom ) )
       ALLOCATE (   vd( MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom ) )
      
       
       Call ReadPDB(PDBid, MinModel, MaxModel, MinChain, MaxChain, MinResidue, MaxResidue, MinAtom, MaxAtom, PDBerror, ATM, AANm, X, Y, Z, B, vd)
      
      
! ------------------------------------------------------------------ Part 0.3 : Read in the PDB Files  [Head Part]
       
       Call HeadPDB(PDBid, PDBRecord, Line)       
       
! ------------------------------------------------------------------ Part 0.5 : Put the coordinates in a sequence
       
       Call SitePDB(PDBid, MinModel, MaxModel, MinChain, MaxChain, MinResidue, MaxResidue, MinAtom, MaxAtom, ATM, AANm, vd, Site)
       
       Allocate ( SeqIndex(1:Site, 1:4) )
       Allocate ( Seq( MinModel:MaxModel, MinChain:MaxChain, MinResidue:MaxResidue, MinAtom:MaxAtom ) )
       
       Call SequPDB(PDBid, MinModel, MaxModel, MinChain, MaxChain, MinResidue, MaxResidue, MinAtom, MaxAtom, ATM, AANm, vd, Site, SeqIndex, Seq)
       
       
!   *   *   *   * 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 *   *   *   *  !       
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
       
       Allocate( xx(1:SITE) )
       Allocate( yy(1:SITE) )
       Allocate( zz(1:SITE) )
       
       Do n=1,SITE
       
           k=SeqIndex(n,1)
           j=SeqIndex(n,2)
           i=SeqIndex(n,3)
           s=SeqIndex(n,4)
           
           xx(n)=X(k,j,i,s)
           yy(n)=Y(k,j,i,s)
           zz(n)=Z(k,j,i,s)
           
       Enddo

        
! Part 4 --------------------------------------------------------------------------- Calculate the Hessian
       Write(9,*) '   '       
       Write(9,'(A,F)') 'The cutoff distance : ', rcut
       Write(9,'(A,F)') 'The force constant  : ', gama
       Write(9,*) '   '
       
       Write(*,*) '   '
       Write(9,*) '   '
       Write(*,'(A)') 'Start generating the Hessian Matrix . . . . . .  ' 
       Write(9,'(A)') 'Start generating the Hessian Matrix . . . . . .  ' 

!   ===========================================================================
!   |         oanm.hes  : for blazpack.x on iGNMtest server                   |
!   | 1xxx_Hessian.dat  : for anmEigen.m as a input file of the Matlab        |
!   ===========================================================================

	   open (40, file=pdbid//'_Hessian.dat')
	   open (41, STATUS='SCRATCH')

       nonzero=0

!      Diagonal clusters of Hessian Matrix : 

       Do j=1,Site !==================================================================== Do j=1,Site
      
!         Initialize Hessian (dn)
          Do s=1,3
             Do t=1,3
               dn(s,t)=0.0
             End Do
          End do

          Do k=1,Site ! +++++++++++++++++++++++++ Do k=1,Site

             bx=xx(j)-xx(k)
             by=yy(j)-yy(k)
             bz=zz(j)-zz(k)
            
             rij = bx*bx + by*by + bz*bz 
            
             If( (j.ne.k) .and. (rij<rcut*rcut) ) then

		        dn(1,1) = dn(1,1) + gama*bx*bx/rij
                dn(2,2) = dn(2,2) + gama*by*by/rij
	            dn(3,3) = dn(3,3) + gama*bz*bz/rij
               
	            dn(1,2) = dn(1,2) + gama*bx*by/rij
	            dn(1,3) = dn(1,3) + gama*bx*bz/rij
	            dn(2,1) = dn(2,1) + gama*by*bx/rij
	            
	            dn(2,3) = dn(2,3) + gama*by*bz/rij
	            dn(3,1) = dn(3,1) + gama*bz*bx/rij
	            dn(3,2) = dn(3,2) + gama*bz*by/rij

             Endif ! If( (j.ne.k) .and. (rij<rcut) )
          
          Enddo ! +++++++++++++++++++++++++++++++++ Do k=1,Site
         
          Do s=1,3
             Do t=1,3
               
               If( ( (3*j+t-3) >= (3*j+s-3) ) .AND. (dn(s,t) .NE. 0.00) ) then 
               
                  nonzero = nonzero + 1
                  
                  write(41,'(2I8,3x,E)') 3*j+s-3, 3*j+t-3, dn(s,t)                  
!                  write(* ,'(2I8,3x,E, 3I8)') 3*j+s-3, 3*j+t-3, dn(s,t) , j, s, t
               Endif
                              
             EndDo ! t=1,3
          EndDo ! s=1,3
          
!          write(41,'(2I8,3x,E)') 3*j, 3*j, dn(3,3)                  
          
	   
      Enddo !==================================================================== Do j=1,Site
      


!     Off-Diagonals clusters of Hessian Matrix : 

       Do j=1,Site !==================================================================== Do j=1,Site

         Do k=1,Site ! +++++++++++++++++++++++++++++++++ Do k=1,Site

            bx = xx(j) - xx(k)
            by = yy(j) - yy(k)
            bz = zz(j) - zz(k)
            
            rij = bx*bx + by*by + bz*bz 
            
            If( (j.ne.k) .and. (rij<rcut*rcut) ) then           

               dn(1,1) = -gama * bx*bx / rij
               dn(2,2) = -gama * by*by / rij
               dn(3,3) = -gama * bz*bz / rij
                                       
               dn(1,2) = -gama * bx*by / rij
               dn(1,3) = -gama * bx*bz / rij
               dn(2,1) = -gama * by*bx / rij
               
               dn(2,3) = -gama * by*bz / rij
               dn(3,1) = -gama * bz*bx / rij
               dn(3,2) = -gama * bz*by / rij

		       Do s=1,3
	              Do t=1,3
	              
                    If( ( (3*k+t-3) >= (3*j+s-3) ) .AND. (dn(s,t) .NE. 0.00) ) then 
                       
                        nonzero = nonzero + 1
                          
                        write(41,'(2I8,3x,E)') 3*j+s-3, 3*k+t-3, dn(s,t)
                                                  
                    Endif
                    
	              Enddo ! t=1,3
		       EndDo ! s=1,3

            EndIf ! If( (j.ne.k) .and. (rij<rcut) ) 
         
         Enddo ! ++++++++++++++++++++++++++++++++++++++++++ Do k=1,Site

       Enddo !============================================================================== Do j=1,Site

       Rewind(41) 
       
       ct=0
       Do k=1, nonzero
         
          Read(41,'(2I8,3x,E)') i, j, hess
         
          If(i<j)  Then
            ct=ct+2
            write(40,'(2I8,3x,E22.15)') i, j, hess
            write(40,'(2I8,3x,E22.15)') j, i, hess
          ElseIf(i==j) Then
            ct=ct+1
            write(40,'(2I8,3x,E22.15)') i, j, hess        
          EndIF
         
       Enddo ! k
       
       nonzero=ct

       Close(40) ! pdbid//'_Hessian.dat'
       Close(41) ! 'oanm.hes'
       
! ------------------------------------------------------------------------------------------------------- Summary
       Spareness = 1.0 * nonzero / (Site*Site*9.0)
       
       Write(*,*) '    '
       Write(9,*) '    '
       Write(*,'(A,I)') 'The total non-zero matrix elements is : ', nonzero
       Write(9,'(A,I)') 'The total non-zero matrix elements is : ', nonzero
       Write(*,*) '    '
       Write(9,*) '    '
       Write(*,'(A,F)') 'The total spareness of the Hessian Matrix is : ', Spareness       
       Write(9,'(A,F)') 'The total spareness of the Hessian Matrix is : ', Spareness       
       Write(*,*) '    '
       Write(9,*) '    '
       
! ------------------------------------------------------------------------------------------------------- Summary
       rijsum = 0 
       ct = 0
       Do i = 1, SITE-1
            Do j = i, SITE
                bx = xx(i) - xx(j)
                by = yy(i) - yy(j)
                bz = zz(i) - zz(j)
                
                rij = bx*bx + by*by + bz*bz 
                
                If(  rij < rcut*rcut ) Then
                    ct = ct + 1
                    rijsum = rijsum + rij
                Endif
            Enddo
       Enddo
       
       avgbond = ct /Float(SITE)
       

       Write(*,'(A,F)') 'The average bond number per node is : ', avgbond 
       Write(*,'(A,F)') 'The average nodes number within the cutoff : ', avgbond*2      
       Write(9,'(A,F)') 'The average bond number per node is : ', avgbond       
       Write(9,'(A,F)') 'The average nodes number within the cutoff : ', avgbond*2      
       Write(9,'(A,F)') 'The sum of rij^2 : ', rijsum      
       Write(*,*) '    '
       Write(9,*) '    '       
       
       Close(9) !  PDBid//'_Hessian.log'        
                     
       INQUIRE ( FILE = 'Network_illustration.inp', EXIST = INPexist )
       If (INPexist==.TRUE.)  Call Network_illustration ( PDBid, rcut, SITE, xx, yy, zz )

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