precompute_periodic2_ewald_data Subroutine

  subroutine precompute_periodic2_ewald_data(plan)
    type(fmm_plan_type), intent(inout) :: plan
    integer(i32) :: img_i, img_j, h1, h2, k_idx
    real(dp) :: alpha, k1, k2, kmag

    call reset_periodic2_ewald_data(plan%periodic_ewald)
    if (.not. use_periodic2_m2l_root_oracle(plan)) return

    alpha = resolve_periodic2_ewald_alpha(plan)
    if (alpha <= 0.0d0) return

    plan%periodic_ewald%axis1 = plan%options%periodic_axes(1)
    plan%periodic_ewald%axis2 = plan%options%periodic_axes(2)
    plan%periodic_ewald%axis_free = 6_i32 - plan%periodic_ewald%axis1 - plan%periodic_ewald%axis2
    if (plan%periodic_ewald%axis_free <= 0_i32 .or. plan%periodic_ewald%axis_free > 3_i32) return

    plan%periodic_ewald%nimg = max(0_i32, plan%options%periodic_image_layers)
    plan%periodic_ewald%kmax = max(1_i32, plan%options%periodic_ewald_layers)
    plan%periodic_ewald%img_outer = plan%periodic_ewald%nimg + plan%periodic_ewald%kmax
    plan%periodic_ewald%alpha = alpha
    plan%periodic_ewald%soft2 = plan%options%softening*plan%options%softening
    plan%periodic_ewald%cell_area = plan%options%periodic_len(1)*plan%options%periodic_len(2)
    if (plan%periodic_ewald%cell_area <= 0.0d0) then
      call reset_periodic2_ewald_data(plan%periodic_ewald)
      return
    end if
    plan%periodic_ewald%k0_pref = 2.0d0*pi_dp/plan%periodic_ewald%cell_area

    plan%periodic_ewald%screen_count = (2_i32*plan%periodic_ewald%img_outer + 1_i32)**2
    plan%periodic_ewald%inner_count = (2_i32*plan%periodic_ewald%nimg + 1_i32)**2
    plan%periodic_ewald%k_count = (2_i32*plan%periodic_ewald%kmax + 1_i32)**2 - 1_i32

    allocate ( &
      plan%periodic_ewald%screen_shift1(plan%periodic_ewald%screen_count), &
      plan%periodic_ewald%screen_shift2(plan%periodic_ewald%screen_count), &
      plan%periodic_ewald%inner_shift1(plan%periodic_ewald%inner_count), &
      plan%periodic_ewald%inner_shift2(plan%periodic_ewald%inner_count) &
      )
    allocate ( &
      plan%periodic_ewald%k1(plan%periodic_ewald%k_count), plan%periodic_ewald%k2(plan%periodic_ewald%k_count), &
      plan%periodic_ewald%kmag(plan%periodic_ewald%k_count), plan%periodic_ewald%karg0(plan%periodic_ewald%k_count), &
      plan%periodic_ewald%kpref1(plan%periodic_ewald%k_count), plan%periodic_ewald%kpref2(plan%periodic_ewald%k_count), &
      plan%periodic_ewald%kprefz(plan%periodic_ewald%k_count) &
      )

    plan%periodic_ewald%screen_count = 0_i32
    do img_i = -plan%periodic_ewald%img_outer, plan%periodic_ewald%img_outer
      do img_j = -plan%periodic_ewald%img_outer, plan%periodic_ewald%img_outer
        plan%periodic_ewald%screen_count = plan%periodic_ewald%screen_count + 1_i32
        plan%periodic_ewald%screen_shift1(plan%periodic_ewald%screen_count) = real(img_i, dp)*plan%options%periodic_len(1)
        plan%periodic_ewald%screen_shift2(plan%periodic_ewald%screen_count) = real(img_j, dp)*plan%options%periodic_len(2)
      end do
    end do

    plan%periodic_ewald%inner_count = 0_i32
    do img_i = -plan%periodic_ewald%nimg, plan%periodic_ewald%nimg
      do img_j = -plan%periodic_ewald%nimg, plan%periodic_ewald%nimg
        plan%periodic_ewald%inner_count = plan%periodic_ewald%inner_count + 1_i32
        plan%periodic_ewald%inner_shift1(plan%periodic_ewald%inner_count) = real(img_i, dp)*plan%options%periodic_len(1)
        plan%periodic_ewald%inner_shift2(plan%periodic_ewald%inner_count) = real(img_j, dp)*plan%options%periodic_len(2)
      end do
    end do

    k_idx = 0_i32
    do h1 = -plan%periodic_ewald%kmax, plan%periodic_ewald%kmax
      k1 = two_pi_dp*real(h1, dp)/plan%options%periodic_len(1)
      do h2 = -plan%periodic_ewald%kmax, plan%periodic_ewald%kmax
        if (h1 == 0_i32 .and. h2 == 0_i32) cycle
        k2 = two_pi_dp*real(h2, dp)/plan%options%periodic_len(2)
        kmag = sqrt(k1*k1 + k2*k2)
        if (kmag <= tiny(1.0d0)) cycle
        k_idx = k_idx + 1_i32
        plan%periodic_ewald%k1(k_idx) = k1
        plan%periodic_ewald%k2(k_idx) = k2
        plan%periodic_ewald%kmag(k_idx) = kmag
        plan%periodic_ewald%karg0(k_idx) = 0.5d0*kmag/alpha
        plan%periodic_ewald%kpref1(k_idx) = pi_dp*k1/(plan%periodic_ewald%cell_area*kmag)
        plan%periodic_ewald%kpref2(k_idx) = pi_dp*k2/(plan%periodic_ewald%cell_area*kmag)
        plan%periodic_ewald%kprefz(k_idx) = pi_dp/plan%periodic_ewald%cell_area
      end do
    end do
    plan%periodic_ewald%k_count = k_idx
    plan%periodic_ewald%ready = .true.
    plan%options%periodic_ewald_alpha = alpha
  end subroutine precompute_periodic2_ewald_data