 program main 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
!
! Simple elastic network model (ENM)
! Protein structure coarse-grained at the amino-acid level
! Residues occupy the alpha-C positions at equilibrium
! All masses equal
! All force constants are equal
! The connectivity graph is built by imposing a cutoff Rc on 
! the equilibrium inter-distances.
!
! The code computes the ENM Hessian, diagonalizes it and 
! prints the NM frequencies and the dusplacement fields of the 
! 1st and last NM modes. 
! The codes used the Lapack routine DSYEV to perform the diaginalization.
!
! Code written for the Class of Scientific computing (3d year physics) 
! at the University of Orleans, France. 
! Instructor: Prof. F. Piazza
! Date released: 11/4/2018
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
!
 INTEGER, PARAMETER :: NX =  2000
 INTEGER, PARAMETER :: LDA = 2000
 INTEGER I, NCA, N, i1, i2, LCA, LWORK, INFO
 
 DOUBLE PRECISION :: xx, yy, zz, dx, dy, dz, RC, u2, Be
 DOUBLE PRECISION, PARAMETER :: PI = 4.d0 * ATAN(1.d0) 
 
 DOUBLE PRECISION, DIMENSION(NX) :: xr,yr,zr    
 DOUBLE PRECISION, DIMENSION(NX) :: BEXP                ! Experimental B factors 
 
 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: WORK  
 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: W     
 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: x,y,z   
 

 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: RR  ! Eq. inter-particle dist. vect
 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: RRm   ! Eq. inter-particle distances
 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: c     ! connectivity matrix
 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: H     ! Hessian matrix

 
 CHARACTER (LEN=100) LINE
 CHARACTER (LEN=8) PDBfile, arg
 CHARACTER (LEN=6) KEY
 CHARACTER (LEN=4) ATOM_TYPE
 
 LOGICAL CONNECT

!==============================================================================
! Read structure 

 if ( iargc() == 2 ) then 
     call getarg(1,PDBfile)
     call getarg(2,arg)
     read (arg,*) Rc     
 else 
     print*, ' Usage: ANM.ex [ PDB_file, Rc ]' 
     STOP
 end if    
 open (1,file=PDBfile)
!==============================================================================
 ic = 0
 NCA = 0
 10 ic = ic + 1 
 read (1,'(A)',end=2000) LINE 
 read (LINE(1:6),'(A6)') KEY
 if ( KEY == 'ATOM  ' ) then
     read(LINE(13:16),'(a4)') ATOM_TYPE  
     read(LINE(31:38),'(f8.3)')  xx 
     read(LINE(39:46),'(f8.3)')  yy
     read(LINE(47:54),'(f8.3)')  zz
     read(LINE(61:66),'(f6.2)')  Be ! Read B-factor    
     if ( ATOM_TYPE == ' CA ') then ! alpha-C found (WARNING: does not handle 
                                    !                structures with multiple 
                                    !                CA/residue )
         NCA = NCA + 1
         if ( NCA > NX ) then 
            print*, ' Currently handles up to ', NX, ' residues'
            STOP
         end if 
         xr(NCA)   = xx; yr(NCA) = yy; zr(NCA) = zz
         BEXP(NCA) = Be
     end if
 end if
 go to 10
 2000 N = ic - 1  
 print*, ' Read coordinates of ', NCA, ' residues'
 
!==============================================================================
! Allocate memory  
 
 LWORK=10*64*9*NCA
 allocate (x(NCA),y(NCA),z(NCA))
 allocate (RR(NCA,NCA,3))
 allocate (RRm(NCA,NCA),c(NCA,NCA))
 allocate (H(LDA,3*NCA)) 
 allocate (W(3*NCA))    
 allocate (WORK(LWORK))
 x = xr(1:NCA); y = yr(1:NCA); z = zr(1:NCA)
!==============================================================================
! Compute inter-particle vectors
 
 RR = 0.d0
 RRm = 0.d0
 do i=1, NCA
    do j=i+1,NCA
       RR(j,i,1) = x(i) - x(j)
       RR(j,i,2) = y(i) - y(j)
       RR(j,i,3) = z(i) - z(j) 
       do ix=1,3     
         RR(i,j,ix) = - RR(j,i,ix)
         RRm(i,j) = RRm(i,j) + RR(i,j,ix)*RR(i,j,ix)
       end do
       RRm(j,i) = RRm(i,j)   ! these are the vector norms squared
       write (20,*) dsqrt(RRm(i,j))
    end do       
 end do
!==============================================================================
! Compute connectivity matrix

 c = 0.d0 
 CONNECT = .false.
 do i=1, NCA
    do j=i+1,NCA
       CONNECT = RRm(i,j) < RC*RC
       if ( CONNECT ) then 
           c(i,j) = 1.d0
           c(j,i) = 1.d0
           CONNECT = .false.
           write (30,'(2i6,2x,f4.1)') i,j
       end if
    end do
 end do
!==============================================================================
! Compute ANM Hessian 

 H = 0.d0
 ! Off-diagonal blocks
 do i=1, NCA
 	do j=i+1, NCA
       do ix=1,3
 	      do iy=1,3
 	      	 i1 = 3*(i-1) + ix
 	      	 i2 = 3*(j-1) + iy
 	      	 H(i1,i2) = - c(i,j)*RR(i,j,ix)*RR(i,j,iy)/RRm(i,j)
 	      	 H(i2,i1) = H(i1,i2) 	     	 
 	      end do
 	   end do    	       	 
 	end do
 end do
 ! Diagonal blocks
 do i=1,NCA
    do ix=1,3
       i1 = 3*(i-1) + ix
       do m=1, NCA
          do iy=1,3	
             i2 = 3*(i-1) + iy
             if ( m .ne. i ) H(i1,i2) = H(i1,i2) +    &
                             c(i,m)*RR(m,i,ix)*RR(m,i,iy)/RRm(m,i)
          end do
       end do      
    end do 
 end do
   
!==============================================================================
! NM analysis. Diagonalize Hessian matrix and write NM frequencies 
!              and selected eigenvectors (NM displacement patterns)

 call DSYEV( 'V', 'U', 3*NCA, H, LDA, W, WORK, LWORK, INFO )
 print*, INFO
 open (1,file='NMfreq.dat')
 do i=1, 3*NCA
    write (1,*) 23.06 * dsqrt(abs(W(i))) ! Frequencies in cm^-1
                                         ! k = 5 kcal/mol/Ang**2
                                         ! M = 110 amu
 end do
 close (1)
 open (1,file='1stHF_NMpattern.dat')     ! Write pattern of 1st NM (HF)
 open (2,file='NthHF_NMpattern.dat')     ! Write pattern of (3N-6)-th NM (LF)
 do i=1, NCA
     write (1,1000) i, (H(3*(i-1)+ix,3*NCA),ix=1,3)
     write (2,1000) i, (H(3*(i-1)+ix,7),ix=1,3)
 end do
 close (1)
 close (2) 
 
!==============================================================================
! Compute mean square displacements
!
!  < u_{ia}^2 > = kB T / m_i \sum_k  \xi^k_{ia} / \omega^2_k
!
!  for particle i along Cartesian direction a
!
!  \xi^k_{ia} : 3N components of k-th eigenvector. These are stored in the k-th 
!               column of the matrix H after call to DSYEV :
!               \xi^k = [ \xi_{1x},\xi_{1y},\xi_{1z}, \xi_{2x},\xi_{2y},\xi_{2z}, 
!                         ....., \xi_{Nx},\xi_{Ny},\xi_{Nz} ]
! \omega^2_k  : The eigenvalues of the dynamical matrix K (REM: masses are all 
!               equal in this simple model) 
!
!
  open (3,file='Temp_factors.dat')     ! Write particles MSDs and exp. B-factors 
  do i=1, NCA
     u2 = 0.d0
     do k=7, 3*NCA
         do ix=1, 3
            u2 = u2 + H(3*(i-1)+ix,k)**2/W(k)
         end do   
     end do
     write (3,'(i5,2(f10.5,4x))') i, 8.d0*PI*u2/3.d0, BEXP(i)
  end do   
  close (3)

!==============================================================================
 1000 format(i5,3f12.6)
!============================================================================== 
 end program



