generate_perturbed_geometries Subroutine

   subroutine generate_perturbed_geometries(reference_geom, displacement, forward_geoms, backward_geoms)
      !! Generate all forward and backward displaced geometries for finite difference calculations
      !!
      !! For a system with N atoms, this generates:
      !!   - 3N forward-displaced geometries (+x, +y, +z for each atom)
      !!   - 3N backward-displaced geometries (-x, -y, -z for each atom)
      !!
      !! These can be used to compute:
      !!   - Gradient: from energies at ±displacement
      !!   - Hessian: from gradients at ±displacement
      !!
      !! Args:
      !!   reference_geom: The reference geometry to perturb
      !!   displacement: Step size in Bohr (typical: 0.005 Bohr)
      !!   forward_geoms: Output array of forward-displaced geometries (size: 3*n_atoms)
      !!   backward_geoms: Output array of backward-displaced geometries (size: 3*n_atoms)
      type(physical_fragment_t), intent(in) :: reference_geom
      real(dp), intent(in) :: displacement
      type(displaced_geometry_t), intent(out), allocatable :: forward_geoms(:)
      type(displaced_geometry_t), intent(out), allocatable :: backward_geoms(:)

      integer :: n_atoms, n_displacements
      integer :: iatom, icoord, idx
      integer :: i

      n_atoms = reference_geom%n_atoms
      n_displacements = 3*n_atoms  ! x, y, z for each atom

      allocate (forward_geoms(n_displacements))
      allocate (backward_geoms(n_displacements))

      ! Generate all displaced geometries
      idx = 0
      do iatom = 1, n_atoms
         do icoord = 1, 3  ! x, y, z
            idx = idx + 1

            ! Forward displacement (+h)
            forward_geoms(idx)%atom_index = iatom
            forward_geoms(idx)%coordinate = icoord
            forward_geoms(idx)%direction = +1
            forward_geoms(idx)%displacement = displacement
            call copy_and_displace_geometry(reference_geom, iatom, icoord, +displacement, &
                                            forward_geoms(idx)%geometry)

            ! Backward displacement (-h)
            backward_geoms(idx)%atom_index = iatom
            backward_geoms(idx)%coordinate = icoord
            backward_geoms(idx)%direction = -1
            backward_geoms(idx)%displacement = displacement
            call copy_and_displace_geometry(reference_geom, iatom, icoord, -displacement, &
                                            backward_geoms(idx)%geometry)
         end do
      end do

   end subroutine generate_perturbed_geometries