OpenMP parallelization for a brownian particle

[AdityaGanesh]

Hi, I am fairly new to OpenMP in fortran. I have the following code where a particle is evolving from the position x = 0 and this is governed by the equation dx = A*zeta*sqrt(dt). where I have two types of zeta, which is a Gaussian random number and a Unit variance uniform random number. As this problem is embarrassingly parallel, I was expecting a scaleup of equal proportions, but however, I couldn't find the code running faster than the optimized (-O3) ifort command. Can you kindly let me know where I can make changes in my code? I have also seen discussions on the random number generator function random_number(x) but was not able to follow it properly. Here is my code. My objective is to increase the value of N:

module parameters

    implicit none

    integer, parameter    :: k = selected_int_kind(9)
    real(8), parameter    :: dt = 0.020d0 ! Time step
    real(8), parameter    :: T = 10.d0    ! Total time
    real(8), parameter    :: A = 1.d0     ! Constant
    integer(k), parameter :: N = 1000000 ! Number of Trajectories

    contains

    !!========================================================================!!
    subroutine Gaussian_random_number(rand)

    !! Generates a Gaussian random number from the uniform random number 
    !! generator using the Box-Muller transform
        implicit none

        real(8), intent(out) :: rand                     ! Subroutine output
        real(8), parameter :: pi = 4.0d0 * atan(1.0d0)   ! Value of pi
        real(8) :: u1, u2                                ! Two random numbers

        !! Generating two uniform random numbers
        call random_number(u1)
        call random_number(u2)

        !! Box-Muller Transform
        rand = sqrt(-2.0d0 * log(u1)) * cos(2.0d0 * pi * u2)

    end subroutine Gaussian_random_number

    !!========================================================================!!
    subroutine Uniform_random_number_with_unit_variance(rand2)
    !! Generates a uniform random number between -sqrt(3) and + sqrt(3)
        implicit none

        real(8), intent(out) :: rand2
        real(8) :: p, q
        p = -sqrt(3.0d0)
        q =  sqrt(3.0d0)
        call random_number(rand2)
        rand2 = p + (q - p) * rand2

    end subroutine Uniform_random_number_with_unit_variance
    !!========================================================================!!

end module parameters

program Brownian_particle
    use parameters
    use omp_lib
    implicit none

    real(8) :: zeta       ! Gaussian random noise
    integer(8) :: Nt      ! Number of time steps
    integer(8) :: i, j    ! Counters
    real(8), allocatable, dimension(:,:) :: x_store
    real(8), allocatable, dimension(:,:) :: x_store_sq
    real(8), allocatable, dimension(:)   :: x_store_sq_avg
    real(8), allocatable, dimension(:)   :: Time
    !$ integer, parameter :: num_threads = 4

    ! For calculating the slope of the given data
    real(8) :: sumx, sumy, sumxx, sumxy, xave, yave, slope

    !$ call omp_set_num_threads(num_threads)

    Nt = T/dt
    allocate(x_store(Nt + 1, N), x_store_sq(Nt + 1, N))
    allocate(x_store_sq_avg(Nt + 1), Time(Nt + 1))
    forall (i = 1 : Nt + 1) Time(i) = (i - 1)*dt
    x_store(1,:) = 0.0d0

    call random_init(.true., .true.)
    !$omp parallel do shared(x_store)
    do j = 1, N
        do i = 2, Nt
            ! call Uniform_random_number_with_unit_variance(zeta) 
            call Gaussian_random_number(zeta)
            x_store(i, j) = x_store(i - 1, j) + A * zeta * sqrt(dt)
        end do
    end do
    !$omp end parallel do
    !$omp barrier

    !c $omp parallel workshare
    x_store_sq = x_store**2
    x_store_sq_avg = sum(x_store_sq, dim=2)/size(x_store_sq, dim=2)
    !c $omp end parallel workshare
    !c omp barrier

    ! open (10, file = 'omp1.dat', status= 'replace')
    ! open(20, file = 'omp2.dat', status='replace')
    !     do i = 0, Nt
    !         write(10, '(1500F14.7)') Time(i+1), (x_store(i, j), j = 1, N)
    !         write(20, '(1000F14.7)')  Time(i+1), x_store_sq_avg(i)
    !     end do
    ! close(10)
    ! close(20)

    ! Calculate the slope using linear regression (y=mx)
    sumx = 0.0d0
    sumy = 0.0d0
    sumxx = 0.0d0
    sumxy = 0.0d0

    !c $omp parallel do reduction(+:sumx, sumy, sumxx, sumxy)
    do i = 0, Nt
        sumx = sumx + x_store_sq_avg(i)
        sumy = sumy + Time(i+1)
        sumxx = sumxx + x_store_sq_avg(i)*x_store_sq_avg(i)
        sumxy = sumxy + x_store_sq_avg(i)*Time(i+1)
    end do
    !c $omp end parallel do
    !c $omp barrier

    xave = sumx / n
    yave = sumy / n
    slope = (sumxy - sumx*yave) / (sumxx - sumx*xave)

    print *, 'The slope of the line is: ', slope

end program Brownian_particle

The commands I am using are as follows

ifort -O3 BP.f90 -o BP -Qopenmp

./BP

I am running this on my windows laptop with oneAPI. My laptop has an AMD Ryzen 5 5500U. The results should print the slope to be ~1.