anchor.diff:              end do
           end do
           riesz_energy = 0.5d0 * riesz_energy
  -    end subroutine riesz_energy
  +    end function riesz_energy
   
   
       pure subroutine riesz_energy_force(points, energy, force)
  
positive.code:                             & image_point, energy, image_force)
                      end do
                  end do
                  grad(:,p) = grad(:,p) + &
                      & matmul(symmetry_group_inv(:,:,a), image_force)
              end do
              grad(:,p) = grad(:,p) - &
                  & dot_product(grad(:,p), points(:,p)) * points(:,p)
          end do
          energy = 0.5d0 * energy
      end subroutine riesz_energy_force
  
negative.code:     pure subroutine add_pair_energy(source_pt, target_pt, energy)
          real(dp), intent(in) :: source_pt(d + 1), target_pt(d + 1)
          real(dp), intent(inout) :: energy
  
          energy = energy + pair_potential(norm2(target_pt - source_pt))
      end subroutine add_pair_energy
  