        !************************************************************************
        !         A Subroutine of Main, Generating the Hessian Matrix           *
        !************************************************************************
       
        Subroutine Hessian_XYZ ( PDBid, rcut, gama )
       
        Implicit None
        
        Integer     SITE       ! Number of vertices
!        Parameter  (SITE=12)   
        
       
! ----- Variables : -----------------------------------------------------------------------!
        Character   PDBid*4
        Real        rcut      ! cutoff distance (in A)
        Real        gama   ! the force constant, gama
        
        
     
       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),    Allocatable  ::   hn(:, :, :, :)  ! Hessian Matrix
       
       Real(8)      avgbond            ! the average bond number per node

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

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

       Open(5, File=PDBid//'.xyz')
       
       read(5,*) SITE
       
       If( SITE<1 ) Then
           Write(*,'(2A)') 'Error! The first line in file [', PDBid//'.xyz] is not the total node number!'
           Write(9,'(2A)') 'Error! The first line in file [', PDBid//'.xyz] is not the total node number!'
           Write(*,*) '   '
           Write(9,*) '   '           
           Goto 999
       EndIf
       
       Write(*,'(3A,I)') 'The total node number read from [', PDBid//'.xyz]', ' file is : ', SITE
       Write(9,'(3A,I)') 'The total node number read from [', PDBid//'.xyz]', ' file is : ', SITE
       Write(*,*) '   '
       Write(9,*) '   '
       
       Allocate( xx(1:SITE) )
       Allocate( yy(1:SITE) )
       Allocate( zz(1:SITE) )
       Allocate( hn(1:SITE, 1:SITE, 1:3, 1:3) ) ! Hessian Matrix
       
       
       ct=0
       Do n=1, SITE
           Read(5, *, ERR=999, END=450) xx(n), yy(n), zz(n)
           ct=ct+1
       Enddo

450    If ( ct==SITE ) Then
           Write(*,'(2A)') 'Successfully read in the coordinates data from [', PDBid//'.xyz] file.'
           Write(9,'(2A)') 'Successfully read in the coordinates data from [', PDBid//'.xyz] file.'
           Write(9,*) '   '
       Else 
           Write(*,'(2A)') 'Error! The actual nodes number in file [', PDBid//'.xyz] is not right!'
           Write(*,'(A,I)') '-- Supposed to read in nodes : ', SITE
           Write(*,'(A,I)') '-- Actually to read in nodes : ', ct
           Write(*,*) '   '
           Write(9,'(2A)') 'Error! The actual nodes number in file [', PDBid//'.xyz] is not right!'
           Write(9,'(A,I)') '-- Supposed to read in nodes : ', SITE
           Write(9,'(A,I)') '-- Actually to read in nodes : ', ct
           Write(9,*) '   '
           Goto 998           
       EndIf

            
       
       Close(5) ! PDBid//'.xyz'


        
! 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        |
!   ===========================================================================

       Do j=1, SITE   ! ------------------- Initialization 
       Do k=1, SITE
           Do s=1,3
           Do t=1,3
               hn(j,k,s,t)=0.0d0
           Enddo
           Enddo
       Enddo
       Enddo


	   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.0d0
             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.0d0) ) 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
          
          
          Do s=1,3
          Do t=1,3
              hn(j,j,s,t) = dn(s,t)
          Enddo
          Enddo
                    
!          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
		       
		       
		       Do s=1,3
		       Do t=1,3
		            hn(j,k,s,t) = dn(s,t) 
		       Enddo
		       Enddo

            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
       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 ) ct = ct + 1
            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(*,*) '    '
       Write(9,*) '    '
       
       
       Close(9) !  PDBid//'_Hessian.log'        
                     
                     
!  +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+
       
!       Open(63, File=PDBID//'_Hessian Matrix.dat')
!        
!        
!        Do j = 1, SITE
!        Do s=1, 3
!        Do k = 1, SITE
!        Do t=1, 3
!            Write(63,*) 3*(j-1)+s, 3*(k-1)+t, hn(j,k,s,t) 
!        Enddo
!        Enddo
!        Enddo ! k
!        Enddo ! j
!        
!       
!       Close(63) ! File=PDBID//'_Hessian Matrix.dat'

998    DeAllocate( xx, yy, zz, hn)

!  +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+ +=+                             

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