src/bob_integrator.f90

! BOB Quantum Civilization Engine - Time Integration
! Module: bob_integrator
! Purpose: Time evolution of quantum states under Hamiltonians
! Standard: Fortran 2018

module bob_integrator
    use bob_kinds
    use bob_errors
    use bob_state
    use bob_hamiltonian
    implicit none
    private
    
    !> Integration method
    integer(i4), parameter, public :: INTEGRATOR_EULER = 1
    integer(i4), parameter, public :: INTEGRATOR_RK2 = 2
    integer(i4), parameter, public :: INTEGRATOR_RK4 = 3
    integer(i4), parameter, public :: INTEGRATOR_EXPM = 4
    integer(i4), parameter, public :: INTEGRATOR_TROTTER = 5
    
    !> Time integrator
    type, public :: bob_time_integrator
        integer(i4) :: method                           ! Integration method
        real(wp) :: dt                                  ! Time step
        real(wp) :: time                                ! Current time
        integer(i8) :: steps_taken                      ! Number of steps
        real(wp) :: error_estimate                      ! Error estimate
        logical(lk) :: adaptive                         ! Adaptive time stepping
        real(wp) :: tolerance                           ! Error tolerance
        character(len=64) :: label                      ! Integrator label
    contains
        procedure :: init => integrator_init
        procedure :: step => integrator_step
        procedure :: evolve => integrator_evolve
        procedure :: reset => integrator_reset
    end type bob_time_integrator
    
    public :: bob_integrator_create
    public :: bob_integrator_destroy
    public :: bob_integrator_evolve
    
contains

    !> Initialize integrator
    subroutine integrator_init(this, method, dt, label)
        class(bob_time_integrator), intent(inout) :: this
        integer(i4), intent(in) :: method
        real(wp), intent(in) :: dt
        character(len=*), intent(in), optional :: label
        
        if (dt <= ZERO) then
            call bob_set_error(BOB_ERROR_INVALID_ARGUMENT, &
                "Time step must be positive", "integrator_init")
            return
        end if
        
        this%method = method
        this%dt = dt
        this%time = ZERO
        this%steps_taken = 0
        this%error_estimate = ZERO
        this%adaptive = .false.
        this%tolerance = 1.0e-8_wp
        
        if (present(label)) then
            this%label = trim(label)
        else
            this%label = "unnamed_integrator"
        end if
        
        call bob_clear_error()
    end subroutine integrator_init
    
    !> Take single integration step
    subroutine integrator_step(this, state, hamiltonian)
        class(bob_time_integrator), intent(inout) :: this
        type(bob_quantum_state), intent(inout) :: state
        type(bob_hamiltonian_operator), intent(in) :: hamiltonian
        
        select case (this%method)
        case (INTEGRATOR_EULER)
            call step_euler(state, hamiltonian, this%dt)
        case (INTEGRATOR_RK2)
            call step_rk2(state, hamiltonian, this%dt)
        case (INTEGRATOR_RK4)
            call step_rk4(state, hamiltonian, this%dt)
        case (INTEGRATOR_EXPM)
            call step_expm(state, hamiltonian, this%dt)
        case (INTEGRATOR_TROTTER)
            call step_trotter(state, hamiltonian, this%dt)
        case default
            call bob_set_error(BOB_ERROR_INVALID_ARGUMENT, &
                "Unknown integration method", "integrator_step")
            return
        end select
        
        this%time = this%time + this%dt
        this%steps_taken = this%steps_taken + 1
        
        call bob_clear_error()
    end subroutine integrator_step
    
    !> Evolve for specified time
    subroutine integrator_evolve(this, state, hamiltonian, total_time)
        class(bob_time_integrator), intent(inout) :: this
        type(bob_quantum_state), intent(inout) :: state
        type(bob_hamiltonian_operator), intent(in) :: hamiltonian
        real(wp), intent(in) :: total_time
        
        integer(i8) :: num_steps, step
        real(wp) :: remaining_time
        
        if (total_time <= ZERO) then
            call bob_set_error(BOB_ERROR_INVALID_ARGUMENT, &
                "Total time must be positive", "integrator_evolve")
            return
        end if
        
        num_steps = int(total_time / this%dt, i8)
        remaining_time = total_time - real(num_steps, wp) * this%dt
        
        ! Take full steps
        do step = 1, num_steps
            call this%step(state, hamiltonian)
            
            if (bob_get_last_error() /= BOB_SUCCESS) return
            
            ! Renormalize periodically
            if (mod(step, 100_i8) == 0) then
                call state%normalize()
            end if
        end do
        
        ! Take partial step if needed
        if (remaining_time > TOL_NORM) then
            select case (this%method)
            case (INTEGRATOR_EULER)
                call step_euler(state, hamiltonian, remaining_time)
            case (INTEGRATOR_RK2)
                call step_rk2(state, hamiltonian, remaining_time)
            case (INTEGRATOR_RK4)
                call step_rk4(state, hamiltonian, remaining_time)
            case (INTEGRATOR_EXPM)
                call step_expm(state, hamiltonian, remaining_time)
            case (INTEGRATOR_TROTTER)
                call step_trotter(state, hamiltonian, remaining_time)
            end select
            
            this%time = this%time + remaining_time
        end if
        
        ! Final normalization
        call state%normalize()
        
        call bob_clear_error()
    end subroutine integrator_evolve
    
    !> Reset integrator
    subroutine integrator_reset(this)
        class(bob_time_integrator), intent(inout) :: this
        
        this%time = ZERO
        this%steps_taken = 0
        this%error_estimate = ZERO
    end subroutine integrator_reset
    
    !> Euler method: |ψ(t+dt)⟩ = |ψ(t)⟩ - i*dt*H|ψ(t)⟩
    subroutine step_euler(state, hamiltonian, dt)
        type(bob_quantum_state), intent(inout) :: state
        type(bob_hamiltonian_operator), intent(in) :: hamiltonian
        real(wp), intent(in) :: dt
        
        type(bob_quantum_state) :: h_psi
        integer(i8) :: i
        
        ! Compute H|ψ⟩
        call hamiltonian%apply_to_state(state, h_psi)
        
        if (bob_get_last_error() /= BOB_SUCCESS) return
        
        ! Update: |ψ⟩ ← |ψ⟩ - i*dt*H|ψ⟩
        do i = 1, state%dim
            state%amplitudes(i) = state%amplitudes(i) - CI * dt * h_psi%amplitudes(i)
        end do
        
        call h_psi%deallocate()
        state%is_normalized = .false.
    end subroutine step_euler
    
    !> Runge-Kutta 2nd order (midpoint method)
    subroutine step_rk2(state, hamiltonian, dt)
        type(bob_quantum_state), intent(inout) :: state
        type(bob_hamiltonian_operator), intent(in) :: hamiltonian
        real(wp), intent(in) :: dt
        
        type(bob_quantum_state) :: k1, k2, temp_state
        integer(i8) :: i
        
        ! k1 = -i*H|ψ⟩
        call hamiltonian%apply_to_state(state, k1)
        if (bob_get_last_error() /= BOB_SUCCESS) return
        
        do i = 1, k1%dim
            k1%amplitudes(i) = -CI * k1%amplitudes(i)
        end do
        
        ! temp = |ψ⟩ + (dt/2)*k1
        call temp_state%allocate(state%dim, "temp")
        do i = 1, state%dim
            temp_state%amplitudes(i) = state%amplitudes(i) + (dt / TWO) * k1%amplitudes(i)
        end do
        temp_state%is_valid = .true.
        
        ! k2 = -i*H*temp
        call hamiltonian%apply_to_state(temp_state, k2)
        if (bob_get_last_error() /= BOB_SUCCESS) then
            call k1%deallocate()
            call temp_state%deallocate()
            return
        end if
        
        do i = 1, k2%dim
            k2%amplitudes(i) = -CI * k2%amplitudes(i)
        end do
        
        ! Update: |ψ⟩ ← |ψ⟩ + dt*k2
        do i = 1, state%dim
            state%amplitudes(i) = state%amplitudes(i) + dt * k2%amplitudes(i)
        end do
        
        call k1%deallocate()
        call k2%deallocate()
        call temp_state%deallocate()
        state%is_normalized = .false.
    end subroutine step_rk2
    
    !> Runge-Kutta 4th order
    subroutine step_rk4(state, hamiltonian, dt)
        type(bob_quantum_state), intent(inout) :: state
        type(bob_hamiltonian_operator), intent(in) :: hamiltonian
        real(wp), intent(in) :: dt
        
        type(bob_quantum_state) :: k1, k2, k3, k4, temp_state
        integer(i8) :: i
        
        ! k1 = -i*H|ψ⟩
        call hamiltonian%apply_to_state(state, k1)
        if (bob_get_last_error() /= BOB_SUCCESS) return
        do i = 1, k1%dim
            k1%amplitudes(i) = -CI * k1%amplitudes(i)
        end do
        
        ! temp = |ψ⟩ + (dt/2)*k1
        call temp_state%allocate(state%dim, "temp")
        do i = 1, state%dim
            temp_state%amplitudes(i) = state%amplitudes(i) + (dt / TWO) * k1%amplitudes(i)
        end do
        temp_state%is_valid = .true.
        
        ! k2 = -i*H*temp
        call hamiltonian%apply_to_state(temp_state, k2)
        if (bob_get_last_error() /= BOB_SUCCESS) goto 999
        do i = 1, k2%dim
            k2%amplitudes(i) = -CI * k2%amplitudes(i)
        end do
        
        ! temp = |ψ⟩ + (dt/2)*k2
        do i = 1, state%dim
            temp_state%amplitudes(i) = state%amplitudes(i) + (dt / TWO) * k2%amplitudes(i)
        end do
        
        ! k3 = -i*H*temp
        call hamiltonian%apply_to_state(temp_state, k3)
        if (bob_get_last_error() /= BOB_SUCCESS) goto 999
        do i = 1, k3%dim
            k3%amplitudes(i) = -CI * k3%amplitudes(i)
        end do
        
        ! temp = |ψ⟩ + dt*k3
        do i = 1, state%dim
            temp_state%amplitudes(i) = state%amplitudes(i) + dt * k3%amplitudes(i)
        end do
        
        ! k4 = -i*H*temp
        call hamiltonian%apply_to_state(temp_state, k4)
        if (bob_get_last_error() /= BOB_SUCCESS) goto 999
        do i = 1, k4%dim
            k4%amplitudes(i) = -CI * k4%amplitudes(i)
        end do
        
        ! Update: |ψ⟩ ← |ψ⟩ + (dt/6)*(k1 + 2*k2 + 2*k3 + k4)
        do i = 1, state%dim
            state%amplitudes(i) = state%amplitudes(i) + &
                (dt / 6.0_wp) * (k1%amplitudes(i) + TWO * k2%amplitudes(i) + &
                                 TWO * k3%amplitudes(i) + k4%amplitudes(i))
        end do
        
999     continue
        call k1%deallocate()
        call k2%deallocate()
        call k3%deallocate()
        call k4%deallocate()
        call temp_state%deallocate()
        state%is_normalized = .false.
    end subroutine step_rk4
    
    !> Matrix exponential method: |ψ(t+dt)⟩ = exp(-i*H*dt)|ψ(t)⟩
    subroutine step_expm(state, hamiltonian, dt)
        type(bob_quantum_state), intent(inout) :: state
        type(bob_hamiltonian_operator), intent(in) :: hamiltonian
        real(wp), intent(in) :: dt
        
        complex(cwp), allocatable :: exp_matrix(:,:)
        complex(cwp), allocatable :: new_amplitudes(:)
        integer(i8) :: i, j, k
        integer :: stat
        real(wp) :: factorial
        complex(cwp) :: term_coeff
        integer, parameter :: MAX_TERMS = 20
        
        ! Allocate matrices
        allocate(exp_matrix(state%dim, state%dim), stat=stat)
        if (stat /= 0) then
            call bob_set_error(BOB_ERROR_ALLOCATION, &
                "Failed to allocate exponential matrix", "step_expm")
            return
        end if
        
        allocate(new_amplitudes(state%dim), stat=stat)
        if (stat /= 0) then
            deallocate(exp_matrix)
            call bob_set_error(BOB_ERROR_ALLOCATION, &
                "Failed to allocate new amplitudes", "step_expm")
            return
        end if
        
        ! Compute exp(-i*H*dt) using Taylor series
        ! exp(A) = I + A + A²/2! + A³/3! + ...
        
        ! Initialize to identity
        exp_matrix = CZERO
        do i = 1, state%dim
            exp_matrix(i,i) = CONE
        end do
        
        ! Add terms
        factorial = ONE
        do k = 1, MAX_TERMS
            factorial = factorial * real(k, wp)
            term_coeff = (-CI * dt) ** k / factorial
            
            ! Add term: (-i*H*dt)^k / k!
            ! Simplified: just add scaled Hamiltonian powers
            do i = 1, state%dim
                do j = 1, state%dim
                    exp_matrix(i,j) = exp_matrix(i,j) + &
                        term_coeff * hamiltonian%matrix(i,j)
                end do
            end do
        end do
        
        ! Apply to state
        new_amplitudes = CZERO
        do i = 1, state%dim
            do j = 1, state%dim
                new_amplitudes(i) = new_amplitudes(i) + &
                    exp_matrix(i,j) * state%amplitudes(j)
            end do
        end do
        
        state%amplitudes = new_amplitudes
        
        deallocate(exp_matrix, new_amplitudes)
        state%is_normalized = .false.
    end subroutine step_expm
    
    !> Trotter decomposition: exp(-i*H*dt) ≈ exp(-i*H₁*dt)exp(-i*H₂*dt)
    subroutine step_trotter(state, hamiltonian, dt)
        type(bob_quantum_state), intent(inout) :: state
        type(bob_hamiltonian_operator), intent(in) :: hamiltonian
        real(wp), intent(in) :: dt
        
        integer(i8) :: i
        complex(cwp) :: phase_factor
        
        ! Simplified Trotter: apply diagonal and off-diagonal parts separately
        
        ! Apply diagonal part: exp(-i*diag(H)*dt)
        do i = 1, state%dim
            phase_factor = exp(-CI * hamiltonian%matrix(i,i) * dt)
            state%amplitudes(i) = state%amplitudes(i) * phase_factor
        end do
        
        ! Off-diagonal part would require more sophisticated treatment
        ! This is a simplified version
        
        state%is_normalized = .false.
    end subroutine step_trotter
    
    !> C ABI: Create integrator
    function bob_integrator_create(method, dt) result(int_ptr) &
        bind(C, name="bob_integrator_create")
        use, intrinsic :: iso_c_binding
        integer(c_int), value :: method
        real(c_double), value :: dt
        type(c_ptr) :: int_ptr
        
        type(bob_time_integrator), pointer :: integrator
        
        allocate(integrator)
        call integrator%init(method, dt)
        int_ptr = c_loc(integrator)
    end function bob_integrator_create
    
    !> C ABI: Destroy integrator
    subroutine bob_integrator_destroy(int_ptr) bind(C, name="bob_integrator_destroy")
        use, intrinsic :: iso_c_binding
        type(c_ptr), value :: int_ptr
        type(bob_time_integrator), pointer :: integrator
        
        if (.not. c_associated(int_ptr)) return
        
        call c_f_pointer(int_ptr, integrator)
        deallocate(integrator)
    end subroutine bob_integrator_destroy
    
    !> C ABI: Evolve state
    function bob_integrator_evolve(int_ptr, state_ptr, ham_ptr, total_time) result(status) &
        bind(C, name="bob_integrator_evolve")
        use, intrinsic :: iso_c_binding
        type(c_ptr), value :: int_ptr, state_ptr, ham_ptr
        real(c_double), value :: total_time
        integer(c_int) :: status
        
        type(bob_time_integrator), pointer :: integrator
        type(bob_quantum_state), pointer :: state
        type(bob_hamiltonian_operator), pointer :: hamiltonian
        
        if (.not. c_associated(int_ptr) .or. &
            .not. c_associated(state_ptr) .or. &
            .not. c_associated(ham_ptr)) then
            status = BOB_ERROR_INVALID_ARGUMENT
            return
        end if
        
        call c_f_pointer(int_ptr, integrator)
        call c_f_pointer(state_ptr, state)
        call c_f_pointer(ham_ptr, hamiltonian)
        
        call integrator%evolve(state, hamiltonian, total_time)
        status = bob_get_last_error()
    end function bob_integrator_evolve

end module bob_integrator

! Made with Bob