anchor.diff:      call print_parameters
       write(*,*)
       call initialize_point_configuration(points, energy, force)
  -#ifndef PCREO_GRAD_DESC
  +#ifdef PCREO_BFGS
       call identity_matrix_4(inv_hess)
   #endif
       write(*,*)
   
       ! TODO: Is there a more natural choice of initial step size?
       step_size = 1.0E-10_rk
  +#ifdef PCREO_TRACK_ANGLE
  +    step_angle = 0.0_rk
  +#endif
   
       iteration_count = 0
       cur_time = current_time()
  
positive.code:     call print_optimization_status
  
      do
  #ifdef PCREO_GRAD_DESC
  
          step_size = quadratic_line_search(points, energy, force, step_size)
          if (step_size == 0.0_rk) then
              call print_optimization_status
              call save_point_file(points, iteration_count)
              write(*,*) "Convergence has been achieved (up to numerical&
                      & round-off error). Exiting."
              stop
          end if
          points = points + step_size * force
          call constrain_points(points)
          call riesz_energy_force(points, energy, new_force)
          step_angle = sum(force * new_force) / &
              & (norm2(force) * norm2(new_force))
          force = new_force
  
  #else
  
  #ifdef PCREO_TRACK_ANGLE
          ! Multiply inverse hessian by force to obtain step direction
  
negative.code:         ! Multiply inverse hessian by negative gradient to obtain step direction
          call matrix_multiply_42(step_direction, inv_hess, force)
  #endif
          step_size = quadratic_line_search( &
                  & points, energy, step_direction, step_size)
          call check_step_size
          new_points = points + step_size * step_direction
          call constrain_points(new_points)
          delta_points = new_points - points
          call riesz_energy_force(new_points, new_energy, new_force)
          delta_gradient = force - new_force
          call bfgs_update_inverse_hessian(inv_hess, delta_points, delta_gradient)
          points = new_points
          energy = new_energy
          force = new_force
  