There is no numerical integration algorithm.
I used double precision so initial error starts from decimal digit 10^(-15).
And final result shows at 10^(1).
Wow.. please be my savior
I used reduction clause and also compared the values between the sum with reduction clause and shared clause.
program 1 - reduction (+:sum)
program 2 - shared (sum0, sum1, sum2) after than sum = sum0 +sum1 +sum2 (0,1 and 2 are thread numbers)
program 1 and 2 show the same result. But different from serial version. (When I use just 1 cpu, it's result is same as serial)
Interestingly, numerical result is different but tendency of particles motion is same.