      program hessians
      implicit none
      character, dimension(20)   :: w
      real, dimension (60)       :: x,y,z,c1,c2,xnew1,ynew1,znew1,xnew2,ynew2,znew2
!xy1, x=derivative with respect to x , y=the coordinates that changes, 1=plus h (2=minus h)
      real, dimension(200)         :: gxx1,gxy1,gxz1,gxx2,gxy2,gxz2,gyx1,gyy1,gyz1,gyx2,gyy2,&
                                      gyz2,gzx1,gzy1,gzz1,gzx2,gzy2,gzz2
      integer                    :: i,j,k,a,b
      double precision,parameter :: h=1.d-5
!f(1,1,i,j,k) x-x component of the second derivative of atom i of diabatical coupling between state j and k
      double precision           :: f(3,3,13,4,4)


      open (unit=60,file="/home/georgia/phenol_yarkony/phenol.4state/geom0",action="read")

      open (unit=70, file="/home/georgia/phenol_yarkony/phenol.4state/x_plus",action="write")
    
      open (unit=80, file="/home/georgia/phenol_yarkony/phenol.4state/x_minus",action="write")
  
      open (unit=90, file="/home/georgia/phenol_yarkony/phenol.4state/y_plus",action="write")

      open (unit=100, file="/home/georgia/phenol_yarkony/phenol.4state/y_minus",action="write")
    
      open (unit=110, file="/home/georgia/phenol_yarkony/phenol.4state/z_plus",action="write")

      open (unit=120, file="/home/georgia/phenol_yarkony/phenol.4state/z_minus",action="write")
   

     10 format (1x,A1,4x,f3.1,2x,4f14.8)
!new geometries to calculate gradients    
      do i = 1,13
       read(60,*,end=100) w(i),c1(i),x(i),y(i),z(i),c2(i)
        xnew1(i)= x(i)+h
        ynew1(i)= y(i)+h
        znew1(i)= z(i)+h
        xnew2(i)= x(i)-h
        ynew2(i)= y(i)-h
        znew2(i)= z(i)-h
       write(70,10)  w(i),c1(i),xnew1(i),y(i),z(i),c2(i)
       write(80,10)  w(i),c1(i),xnew2(i),y(i),z(i),c2(i)
       write(90,10)  w(i),c1(i),x(i),ynew1(i),z(i),c2(i)
       write(100,10) w(i),c1(i),x(i),ynew2(i),z(i),c2(i)
       write(110,10) w(i),c1(i),x(i),y(i),znew1(i),c2(i)
       write(120,10) w(i),c1(i),x(i),y(i),znew2(i),c2(i)
      end do
      
      100 close(60)
          close(70)
          close(80)
          close(90)
          close(100)
          close(110)
          close(120)


      call  system("cd ~/phenol_yarkony/phenol.4state; ../evalsurf/source/evalsurf.x x_plus >output_x_plus")
      call  system("cd ~/phenol_yarkony/phenol.4state; ../evalsurf/source/evalsurf.x x_minus >output_x_minus")
      call  system("cd ~/phenol_yarkony/phenol.4state; ../evalsurf/source/evalsurf.x y_plus >output_y_plus")
      call  system("cd ~/phenol_yarkony/phenol.4state; ../evalsurf/source/evalsurf.x y_minus >output_y_minus")
      call  system("cd ~/phenol_yarkony/phenol.4state; ../evalsurf/source/evalsurf.x z_plus >output_z_plus")
      call  system("cd ~/phenol_yarkony/phenol.4state; ../evalsurf/source/evalsurf.x z_minus >output_z_minus")

      open (unit=130,file="/home/georgia/phenol_yarkony/phenol.4state/output_x_plus",action="read")
      open (unit=140,file="/home/georgia/phenol_yarkony/phenol.4state/output_x_minus",action="read")
      open (unit=150,file="/home/georgia/phenol_yarkony/phenol.4state/output_y_plus",action="read")
      open (unit=160,file="/home/georgia/phenol_yarkony/phenol.4state/output_y_minus",action="read")
      open (unit=170,file="/home/georgia/phenol_yarkony/phenol.4state/output_z_plus",action="read")
      open (unit=180,file="/home/georgia/phenol_yarkony/phenol.4state/output_z_minus",action="read")

!just reading the output file of the spefic program until we reach the gradients results     
       do i=1,253
          read(130,*)
           read(140,*)
            read(150,*)
             read(160,*)
              read(170,*)
               read(180,*)
        end do
  
          do k=1,10  ! 10 combinations of coupling between states e.g. 11,12,13,43
             do j=1,17 !skip the adiabatic gradients of yarkony program
              read(130,*)
               read(140,*)
                read(150,*)
                 read(160,*)
                  read(170,*)
                   read(180,*)
               end do
          
          do i=1,13
             read(130,*) gxx1(i),gxy1(i),gxz1(i)
               read(140,*) gxx2(i),gxy2(i),gxz2(i)
                 read(150,*) gyx1(i),gyy1(i),gyz1(i)
                   read(160,*) gyx2(i),gyy2(i),gyz2(i)
                     read(170,*) gzx1(i),gzy1(i),gzz1(i)
                        read(180,*) gzx2(i),gzy2(i),gzz2(i)

           if (k.le.4) then 
             a=1
             b=k
          else if (k.le.7) then
            a=2
            b=k-3
          else if (k.eq.9) then
             a=3
             b=k-5
          else 
             a=4
             b=4
          endif

!calculate the hessian matrix
              f(1,1,i,a,b) = (gxx1(i)-gxx2(i))/2*h
              f(2,2,i,a,b) = (gyy1(i)-gyy2(i))/2*h
              f(3,3,i,a,b) = (gzz1(i)-gzz2(i))/2*h
              f(1,2,i,a,b) = (gyx1(i)-gyx2(i))/2*h
              f(2,1,i,a,b) = f(1,2,i,a,b)
              f(1,3,i,a,b) = (gzx1(i)-gzx2(i))/2*h
              f(3,1,i,a,b) = f(1,3,i,a,b)
              f(2,3,i,a,b) = (gzy1(i)-gzy2(i))/2*h
              f(3,2,i,a,b) = f(2,3,i,a,b)
                    
        enddo
    enddo
    
    close (130)
    close (140)
    close (150)
    close (160)
    close (170)
    close (180)

end program hessians
