xtb_calc_gradient Subroutine

   subroutine xtb_calc_gradient(this, fragment, result)
      !! Calculate energy gradient using Extended Tight-Binding (xTB) method
      class(xtb_method_t), intent(in) :: this
      type(physical_fragment_t), intent(in) :: fragment
      type(calculation_result_t), intent(out) :: result

      ! tblite calculation variables
      type(error_type), allocatable :: error
      type(structure_type) :: mol
      real(wp), allocatable :: xyz(:, :)
      integer, allocatable :: num(:)
      type(xtb_calculator) :: calc
      type(wavefunction_type) :: wfn
      real(wp) :: energy
      type(context_type) :: ctx
      integer :: verbosity
      real(wp), allocatable :: gradient(:, :)
      real(wp), allocatable :: sigma(:, :)
      real(wp) :: dipole_wp(3)

      if (this%verbose) then
         print *, "XTB: Calculating gradient using ", this%variant
         print *, "XTB: Fragment has", fragment%n_atoms, "atoms"
         print *, "XTB: nelec =", fragment%nelec
         print *, "XTB: charge =", fragment%charge
         if (allocated(this%solvent)) then
            if (allocated(this%solvation_model)) then
               print *, "XTB: Solvation: ", trim(this%solvation_model), " with solvent = ", this%solvent
            else
               print *, "XTB: Solvation: alpb with solvent = ", this%solvent
            end if
         else
            print *, "XTB: Solvation: none (gas phase)"
         end if
      end if

      ! Convert fragment to tblite format
      allocate (num(fragment%n_atoms))
      allocate (xyz(3, fragment%n_atoms))

      num = fragment%element_numbers
      xyz = fragment%coordinates  ! Already in Bohr

      ! Create molecular structure
      call new(mol, num, xyz, charge=real(fragment%charge, wp), &
               uhf=fragment%multiplicity - 1)

      ! Select and create appropriate GFN calculator
      select case (this%variant)
      case ("gfn1")
         call new_gfn1_calculator(calc, mol, error)
      case ("gfn2")
         call new_gfn2_calculator(calc, mol, error)
      case default
         call result%error%set(ERROR_VALIDATION, "Unknown XTB variant: "//this%variant)
         result%has_error = .true.
         return
      end select

      if (allocated(error)) then
         call result%error%set(ERROR_GENERIC, "Failed to create XTB calculator")
         result%has_error = .true.
         return
      end if

      ! Add solvation if configured (either solvent name or direct dielectric)
      if (allocated(this%solvent) .or. this%dielectric > 0.0_wp) then
         if (allocated(this%solvation_model)) then
            call add_solvation_to_calc(calc, mol, this%solvent, this%solvation_model, this%variant, &
                                       this%use_cds, this%use_shift, this%dielectric, &
                                       this%cpcm_nang, this%cpcm_rscale, error)
         else
            call add_solvation_to_calc(calc, mol, this%solvent, "alpb", this%variant, &
                                       this%use_cds, this%use_shift, this%dielectric, &
                                       this%cpcm_nang, this%cpcm_rscale, error)
         end if
         if (allocated(error)) then
            call result%error%set(ERROR_GENERIC, "Failed to add solvation: "//error%message)
            result%has_error = .true.
            return
         end if
      end if

      ! Allocate gradient and sigma arrays (initialize to zero)
      allocate (gradient(3, fragment%n_atoms))
      allocate (sigma(3, 3))
      gradient = 0.0_wp
      sigma = 0.0_wp

      ! Create wavefunction and run single point calculation with gradient
      call new_wavefunction(wfn, mol%nat, calc%bas%nsh, calc%bas%nao, 1, this%kt, grad=.true.)
      energy = 0.0_wp

      verbosity = merge(1, 0, this%verbose)
      call xtb_singlepoint(ctx, mol, calc, wfn, this%accuracy, energy, &
                           gradient=gradient, sigma=sigma, verbosity=verbosity)

      ! Compute molecular dipole moment from wavefunction
      dipole_wp(:) = matmul(mol%xyz, wfn%qat(:, 1)) + sum(wfn%dpat(:, :, 1), 2)

      ! Store results (XTB is a semi-empirical method, store as SCF energy)
      result%energy%scf = real(energy, dp)
      result%has_energy = .true.

      ! Store gradient
      allocate (result%gradient(3, fragment%n_atoms))
      result%gradient = real(gradient, dp)
      result%has_gradient = .true.

      ! Store sigma (stress tensor)
      allocate (result%sigma(3, 3))
      result%sigma = real(sigma, dp)
      result%has_sigma = .true.

      ! Store dipole moment
      allocate (result%dipole(3))
      result%dipole = real(dipole_wp, dp)
      result%has_dipole = .true.

      if (this%verbose) then
         print *, "XTB: Energy =", result%energy%total()
         print *, "XTB: Gradient norm =", sqrt(sum(result%gradient**2))
         print *, "XTB: Dipole (e*Bohr) =", result%dipole
         print *, "XTB: Dipole magnitude (Debye) =", norm2(result%dipole)*2.541746_dp
         print *, "XTB: Gradient calculation complete"
      end if

      deallocate (num, xyz, gradient, sigma)

   end subroutine xtb_calc_gradient