What are the requirements to do a OpenMP reduction in Fortran?

Question

I am doing molecular dynamics simulation, and trying to parallelize my serial code with the help of openmp constructs in fortran. At, first I am just trying to parallelize the force subroutine, which has a nested do loop in my program. Now I have successfully calculated the i'th and j'th force but when I am using reduction clause to reduce the value of force on an array for i'th molecule over all j'th molecule then it is showing some wrong results. Though my whole molecular dynamics code is too long. But I am providing the force subroutine below.

In this code, when calculating the potential, then every vij between i'th and j'th molecule is calculating in a correct way, then when it is being added to vt while using reduction clausew it gives correct value for total potential, but when calculating fxij is fine,but reduced flx array is not.

subroutine omp_force

        do  i = 1,n    
          flx(i,1) = 0.0d0          
        enddo
    
    !$omp parallel shared(stepmd,vt,flx) 
    !$omp do private(rij,rxij,ryij,rzij,rxijsq,ryijsq,rzijsq,rijsq) &
    !$omp private(xmn,exmn,rng,erng,epsnew,vij) &
    !$omp private(rho,dudr,dudphi,dudphie,fconst,fconste,kappadotr,kappadotre) &
    !$omp private(rcx1,ercx1,rcy1,ercy1,rcz1,ercz1,f_1x,f_1y,f_1z,f_2x,f_2y,f_2z,f_3x,f_3y,f_3z) &
    !$omp private(fxij,fyij,fzij) &
    !$omp schedule(dynamic) reduction(+:vt,flx) 
    
    do i=1,n-1
         ia=i
      do j=i+1,n
           jb=j
           
        rho   = (rij - rng + sigmin)/sigmin
        vij = (1.0d0/rho)**12-(1.0d0/rho)**6    ! calculation of potential between i'th and j'th mol
        
        vt = vt + 4.0d0*(vij)*epsnew
        
        ! total force
        
        fxij = - f_1x + f_2x - f_3x ! force b/w i'th and j'th molecule  

        flx(i,1) = flx(i,1) + fxij
        flx(j,1) = flx(j,1) - fxij
        flx(j,1) = flx(j,1) - fxij
                      end if
                  end if
              end if
          end do
        end do
        
        !$omp end do
        !$omp end parallel
        
        do i=1,n
        write(*,*)'calculates the force',flx(i,1),'for',i !'thread',omp_get_thread_num(),
        end do
end subroutine omp_force