      program yarkony
      implicit none
      character, dimension(20)     :: w
      real, dimension (60)         :: c1,c2,x1,y1,z1,x2,y2,z2,x,y,z
      real, dimension(200)         :: gx1,gy1,gz1,gx2,gy2,gz2
      integer                      :: i,j,k,a,b,c,m,l,ii,a1,a2,a3,ilbl,jlbl,klbl,k1lbl,k2lbl,k3lbl,k4lbl,k5lbl,kk,fl,&
      k6lbl,k7lbl,k8lbl,k9lbl,k10lbl,k11lbl,k14lbl,k15lbl,k16lbl,k17lbl,k18lbl,k19lbl
      double precision,parameter   :: h=1.d-5
      real, dimension(3,13,4,4)    :: deriv1
      real, dimension(39,39,4,4)   :: deriv2
      real, dimension(4,4)         :: u   
      real, dimension(4)           :: v
      real, dimension(39)          :: gpoint
      real, dimension (13,3)       :: coo
      character(len=200)           :: filename,filepath,exl1,exl2,exl3,filename1,filename2,filename3,exl4,exl5,&
      filename4,filename5,command,prog,filepath1,filepath2,filename6,filename7,filename8,filename9,filename10,&
      filepath6,filepath7,filepath9,filepath10,filepath3,exl6,exl7,exl8,exl9,exl10,exl11
      real                         :: dv

      20 format (1x,A1,4x,f3.1,2x,4f14.8)
      
      CALL GET_COMMAND_ARGUMENT(1, FileName)
      fl=index(filename,' ')
     open (unit=10,file=filename, status="unknown",access="sequential", form="formatted", action="read")   
        do i=1,39
          read(10,*,end=100) gpoint(i)
        enddo
      100 close(10)
      a1=1
      a2=2
      a3=3
       do i = 1,13
        x(i)=gpoint(a1)
        y(i)=gpoint(a2)
        z(i)=gpoint(a3)
        a1=a1+3
        a2=a2+3
        a3=a3+3
       enddo
 
    open (unit=10,file=filename, status="unknown",action="write")

    filepath = "/home/gchri/phenol/yarkony/yarkony_prog/phenol.4state"
    ilbl=index(filepath,' ')
    prog = "../evalsurf/source/evalsurf.x"
    jlbl=index(prog,' ')
    command = "cd ; "
   
    filename1=filename(1:fl-1)//"_geom0"
    k1lbl=index(filename1,' ')
    filepath1 = filepath(1:ilbl-1)//"/"//filename1(1:k1lbl)
     k11lbl=index(filepath1,' ')
    exl4= "cp "//filepath(1:ilbl-1)//"/geom0 "//filepath1(1:k11lbl)
    call execute_command_line(exl4) 

    open (unit=60, file=filepath1, status="unknown", action="read")
 
    filename2=filename(1:fl-1)//"_geom_carte"
    k2lbl=index(filename2,' ')
    filepath2 = filepath(1:ilbl-1)//"/"//filename2(1:k2lbl)

    open (unit=40,file=filepath2, status="unknown", action="write")

       do i = 1,13
        read(60,*)  w(i),c1(i),coo(i,1),coo(i,2),coo(i,3),c2(i)
        write(40,20) w(i),c1(i),x(i),y(i),z(i),c2(i) 
       enddo
          close(60)
          close(40)
  
    filename8=filename(1:fl-1)//"_output_dev1"
    k8lbl=index(filename8,' ',.true.)
    filename3 =filename2(1:k2lbl)// " > "//filename8(1:k8lbl)
    k3lbl=index(filename3,' ',.true.)
    exl1 = command(1:3)//filepath(1:ilbl-1)//command(4:5)//prog(1:jlbl+1)//filename3(1:k3lbl)
    filepath3=filepath(1:ilbl-1)//"/"//filename8(1:k8lbl)
    call execute_command_line(exl1)

    open (unit=50,file=filepath3, status="old",action="read")


        do i=1,247
         read(50,*,end=500)
        end do
 
        read(50,*,end=500)  u(1,1), u(1,2), u(1,3), u(1,4)
        read(50,*,end=500) u(2,1), u(2,2), u(2,3), u(2,4)
        read(50,*,end=500) u(3,1), u(3,2), u(3,3), u(3,4)
        read(50,*,end=500) u(4,1), u(4,2), u(4,3), u(4,4)
      
        read(50,*) 
        read(50,*) v(1), v(2), v(3), v(4)
       
       write(10,*) v(1), u(1,2), u(1,3), u(1,4)
       write(10,*) u(2,1), v(2), u(2,3), u(2,4)
       write(10,*) u(3,1), u(3,2), v(3), u(3,4)
       write(10,*) u(4,1), u(4,2), u(4,3), v(4)
      
         do i=1,2
         read(50,*,end=500)
        end do

       do k=1,10
           if (k.le.4) then
             a=1
             b=k
          else if (k.le.7) then
            a=2
            b=k-3
          else if (k.le.9) then
             a=3
             b=k-5
          else
             a=4
             b=4
          endif
            dv=v(b)-v(a)
          if (a.eq.b) then
           do i=1,13
            read(50,*,end=500) deriv1(1,i,a,b),deriv1(2,i,a,b),deriv1(3,i,a,b)
            write(10,*) deriv1(1,i,a,b),deriv1(2,i,a,b),deriv1(3,i,a,b)
           enddo
          else
            do i=1,13
             read(50,*,end=500) deriv1(1,i,a,b),deriv1(2,i,a,b),deriv1(3,i,a,b)
             deriv1(1,i,a,b)=deriv1(1,i,a,b)*dv
             deriv1(2,i,a,b)=deriv1(2,i,a,b)*dv
             deriv1(3,i,a,b)=deriv1(3,i,a,b)*dv
             write(10,*) deriv1(1,i,a,b),deriv1(2,i,a,b),deriv1(3,i,a,b)
            enddo
          endif
            do j=1,17 !skip the diabatic gradients of yarkony program
               read(50,*,end=500)
              enddo
       enddo
      500 close(50)


        do k=1,10
           if (k.le.4) then
             a=1

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

     filename6=filename(1:fl-1)//"_geom_plus"
     k6lbl=index(filename6,' ')
     filepath6 = filepath(1:ilbl-1)//"/"//filename6(1:k6lbl)
    open (unit=70, file=filepath6, status="unknown", action="write")

     filename7=filename(1:fl-1)//"_geom_minus"
    k7lbl=index(filename7,' ')
    filepath7 = filepath(1:ilbl-1)//"/"//filename7(1:k7lbl)
  open (unit=80, file=filepath7, status="unknown", action="write")

             if (i.le.13) then
               do l=1,13
                 if (l.eq.i) then
                   x1(l)=x(l)+h
                   x2(l)=x(l)-h
                 else
                  x1(l)=x(l)
                  x2(l)=x(l)
                 endif
                  y1(l)=y(l)
                  y2(l)=y(l)
                  z1(l)=z(l)
                  z2(l)=z(l)
                write(70,20)  w(l),c1(l),x1(l),y1(l),z1(l),c2(l)
                write(80,20)  w(l),c1(l),x2(l),y2(l),z2(l),c2(l)
               enddo

              else if (i.le.26) then 
                  ii=i-13
                do l=1,13
                  if (l.eq.ii) then
                  y1(l)=y(l)+h
                  y2(l)=y(l)-h
                 else
                  y1(l)=y(l)
                  y2(l)=y(l)
                endif
                  x1(l)=x(l)
                  x2(l)=x(l)
                  z1(l)=z(l)
                  z2(l)=z(l)
               write(70,20)  w(l),c1(l),x1(l),y1(l),z1(l),c2(l)
               write(80,20)  w(l),c1(l),x2(l),y2(l),z2(l),c2(l)
              enddo
            else 
                  ii=i-26
                do l=1,13
                  if (l.eq.ii) then
                  z1(l)=z(l)+h
                  z2(l)=z(l)-h
                else
                  z1(l)=z(l)
                  z2(l)=z(l)
                endif
                  x1(l)=x(l)
                  x2(l)=x(l)
                  y1(l)=y(l)
                  y2(l)=y(l)
               write(70,20)  w(l),c1(l),x1(l),y1(l),z1(l),c2(l)
               write(80,20)  w(l),c1(l),x2(l),y2(l),z2(l),c2(l)
             enddo
           endif
    
    filename9=filename(1:fl-1)//"_output_plus"
    k9lbl=index(filename9,' ',.true.)
    filename4 = filename6(1:k6lbl)//" > "//filename9(1:k9lbl)
    k4lbl=index(filename4,' ',.true.)
    exl2 = command(1:3)//filepath(1:ilbl-1)//command(4:5)//prog(1:jlbl+1)//filename4(1:k4lbl)
      
    filename10=filename(1:fl-1)//"_output_minus"
    k10lbl=index(filename10,' ',.true.)
    filename5 = filename7(1:k7lbl)//" > "//filename10(1:k10lbl)
    k5lbl=index(filename5,' ',.true.)
    exl3 = command(1:3)//filepath(1:ilbl-1)//command(4:5)//prog(1:jlbl+1)//filename5(1:k4lbl)
 
    call execute_command_line(exl2)
    call execute_command_line(exl3)

    filepath9 = filepath(1:ilbl-1)//"/"//filename9(1:k9lbl)
    filepath10 = filepath(1:ilbl-1)//"/"//filename10(1:k10lbl)
   
     open (unit=90,file=filepath9, status="old", action="read")
     open (unit=100,file=filepath10, status="old", action="read")
    

    kk=255+((k-1)*30)
     do m=1,kk
         read(90,*,end=200)
         read(100,*,end=300)
      enddo

        do l=1,13
           read(90,*,end=200) gx1(l),gy1(l),gz1(l)
           read(100,*,end=300) gx2(l),gy2(l),gz2(l)
         enddo
          do m=1,39
            dv=v(b)-v(a)
            if (a.eq.b) then
               if (m.le.13) then
                deriv2(i,m,a,b) = (gx1(m)-gx2(m))/(2*h)
              else if (m.le.26) then
               ii=m-13
               deriv2(i,m,a,b) = (gy1(ii)-gy2(ii))/(2*h)
              else
               ii=m-26
               deriv2(i,m,a,b) = (gz1(ii)-gz2(ii))/(2*h)
             endif 
            else
             if (m.le.13) then
                deriv2(i,m,a,b) = ((gx1(m)-gx2(m)*dv))/(2*h)
              else if (m.le.26) then
               ii=m-13
               deriv2(i,m,a,b) = ((gy1(ii)-gy2(ii))*dv)/(2*h)
              else
               ii=m-26
               deriv2(i,m,a,b) = ((gz1(ii)-gz2(ii))*dv)/(2*h)
             endif
            endif 
          enddo
     close (70)
     close (80)
     200 close (90)
     300 close (100)
   enddo
 enddo
       
 do a=1,4
    do b=a,4
      do i=1,39
        write(10,*)(deriv2(i,j,a,b), j=1,39)
       enddo
    enddo
 enddo
  close(10)

    k14lbl=index(filepath2,' ',.true.)
    k15lbl=index(filepath3,' ',.true.)
    k16lbl=index(filepath6,' ',.true.)
    k17lbl=index(filepath7,' ',.true.)
    k18lbl=index(filepath9,' ',.true.)
    k19lbl=index(filepath10,' ',.true.)   

    exl5 = "rm "//filepath1(1:k11lbl)
    exl6 = "rm "//filepath2(1:k14lbl)
    exl7 = "rm "//filepath3(1:k15lbl)
    exl8 = "rm "//filepath6(1:k16lbl)
    exl9 = "rm "//filepath7(1:k17lbl)
    exl10 = "rm "//filepath9(1:k18lbl)
    exl11 = "rm "//filepath10(1:k19lbl)
    
    call execute_command_line(exl5)
    call execute_command_line(exl6)
    call execute_command_line(exl7)
    call execute_command_line(exl8)
    call execute_command_line(exl9)
    call execute_command_line(exl10)
    call execute_command_line(exl11)




  end program yarkony
