turbomole_gradient Subroutine

subroutine turbomole_gradient(mol, fname, energy, gradient, stat)
   type(structure_type),intent(in) :: mol
   character(len=*),intent(in) :: fname
   real(wp),intent(in) :: energy
   real(wp),intent(in) :: gradient(:, :)
   integer, intent(out) :: stat
   character(len=:),allocatable :: line
   integer  :: i,icycle,line_number
   integer  :: err
   integer  :: igrad ! file handle
   logical  :: exist
   real(wp) :: escf
   real(wp),allocatable :: gscf(:,:)
   real(wp),allocatable :: xyz (:,:)
   allocate( gscf(3,mol%nat), source = 0.0_wp )
   stat = 0
   icycle = 1
   i = 0
   escf = 0.0_wp
   line_number = 0

   inquire(file=fname,exist=exist)
   if (exist) then
      open(newunit=igrad,file=fname)
      read_file: do
         call getline(igrad,line,iostat=err)
         if (err/=0) exit read_file
         i=i+1
         if (index(line,"cycle") > 0) line_number = i
      end do read_file
      if (line_number < 2) then
         stat = 1
         return
      end if

      rewind(igrad)
      skip_lines: do i = 1, line_number-1
         read(igrad,"(a)")
      end do skip_lines
      call getline(igrad,line)
      read(line(10:17),*,iostat=err) icycle
      if (err == 0) read(line(33:51),*,iostat=err) escf
      if (err /= 0) then
         stat = 1
         return
      end if

      allocate(xyz(3,mol%nat))
      do i = 1, mol%nat
         call getline(igrad,line)
         read(line,*,iostat=err) xyz(1,i),xyz(2,i),xyz(3,i)
         if (err /= 0) then
            stat = 1
            return
         end if
      end do
      if (any(abs(xyz-mol%xyz) > 1.0e-8_wp)) then
         stat = 1
         return
      end if
      do i = 1, mol%nat
         call getline(igrad,line)
         read(line,*,iostat=err) gscf(1,i),gscf(2,i),gscf(3,i)
         if (err /= 0) then
            stat = 1
            return
         end if
      end do
      do i = 1, mol%nat
         backspace(igrad)
         backspace(igrad)
      end do
      backspace(igrad)
   else
      open(newunit=igrad,file=fname)
      write(igrad,'("$grad")')
   end if

   write(igrad,'(2x,"cycle =",1x,i6,4x,"SCF energy =",f18.11,3x,'//&
                   '"|dE/dxyz| =",f10.6)') &
      icycle, energy+escf, norm2(gradient+gscf)
   do i = 1, mol%nat
      write(igrad,"(3(F20.14,2x),4x,a2)") mol%xyz(1,i),mol%xyz(2,i),mol%xyz(3,i),mol%sym(i)
   end do
   do i = 1, mol%nat
      write(igrad,"(3D22.13)") gradient(1,i)+gscf(1,i),gradient(2,i)+gscf(2,i),gradient(3,i)+gscf(3,i)
   end do
   write(igrad,'("$end")')
   close(igrad)

end subroutine turbomole_gradient