Сводка OpenMP и массива с Fortran 90

Question

Я пытаюсь вычислить тензор давления кристаллической структуры. Чтобы сделать это, я должен идти Повсеместно все пары частиц, как в упрощать код ниже

do i=1, atom_number  ! sum over atoms i 
    type1 = ATOMS(i) 
    do nj=POINT(i), POINT(i+1)-1 ! sum over atoms j of i's atoms list 
     j = LIST(nj) 
     type2 = ATOMS(j) 
     call get_ff_param(type1,type2,Aab,Bab,Cab,Dab) 
     call distance_sqr2(i,j,r,VECT_R) 
     call gettensor_rij(i,j,T) 
     r = sqrt(r) 
     if (r<coub_cutoff) then 
      local_sum_real(id+1) = local_sum_real(id+1) + (erfc(alpha_ewald*r) 
      derivative = -(erfc(alpha_ewald*r) 
      ff = - (1.0d0/r)*derivative 
      STRESS_EWALDD = STRESS_EWALDD + ff * T ! A 3x3 matrix  
      FDX(i) = FDX(i) + VECT_R(1) * ff 
      FDY(i) = FDY(i) + VECT_R(2) * ff 
      FDZ(i) = FDZ(i) + VECT_R(3) * ff 
      FDX(j) = FDX(j) - VECT_R(1) * ff 
      FDY(j) = FDY(j) - VECT_R(2) * ff 
      FDZ(j) = FDZ(j) - VECT_R(3) * ff  
     end if         
    end do    ! sum over atoms j 

    sum_kin = sum_kin + ATMMASS(i) * (VX(i)**2 + VY(i)**2 + VZ(i)**2) 
    call gettensor_v(i,Q) 
    STRESS_KINETIC = STRESS_KINETIC + ATMMASS(i) * Q 

end do 

Я пытался распараллеливание этой двойной петли с директивой «paralle делать», но есть некоторые проблемы с тензором STRESS_EWALDD и силы FDX .... Поэтому я попытался назначить вручную каждую нить нескольким частицам, как в коде, но все же получить неправильное значение тензора.

!$omp parallel private(id,j,i,nj,type1,type2,T,Q,r,VECT_R,derivative,Aab,Bab,Cab,Dab,ff) 

id = omp_get_thread_num()  
do i=id+1, atom_number,nthreads  ! sum over atoms i 
     type1 = ATOMS(i) 
     do nj=POINT(i), POINT(i+1)-1 ! sum over atoms j of i's atoms list 
      j = LIST(nj) 
      type2 = ATOMS(j) 
      call get_ff_param(type1,type2,Aab,Bab,Cab,Dab) 
      call distance_sqr2(i,j,r,VECT_R) 
      call gettensor_rij(i,j,T) 
      r = sqrt(r) 
      if (r<coub_cutoff) then 
       local_sum_real(id+1) = local_sum_real(id+1) + (erfc(alpha_ewald*r) 
       derivative = -(erfc(alpha_ewald*r) 
       ff = - (1.0d0/r)*derivative 
      local_STRESS_EWALDD(id+1,:,:) = local_STRESS_EWALDD(id+1,:,:) + ff * T  
      FDX(i) = FDX(i) + VECT_R(1) * ff 
       FDY(i) = FDY(i) + VECT_R(2) * ff 
       FDZ(i) = FDZ(i) + VECT_R(3) * ff 
       FDX(j) = FDX(j) - VECT_R(1) * ff 
       FDY(j) = FDY(j) - VECT_R(2) * ff 
       FDZ(j) = FDZ(j) - VECT_R(3) * ff  
      end if         
     end do    ! sum over atoms j 

    local_sum_kin(id+1) = local_sum_kin(id+1) + ATMMASS(i) * (VX(i)**2 + VY(i)**2 + VZ(i)**2) 
    call gettensor_v(i,Q) 
    local_STRESS_KINETIC(id+1,:,:) = local_STRESS_KINETIC(id+1,:,:) + ATMMASS(i) * Q 

end do    ! sum over atoms i 

!$omp end parallel 

do i=1,nthreads 
    sum_real = sum_real + local_sum_real(i) 
    sum_virial_buck = sum_virial_buck + local_sum_virial_buck(i) 
    STRESS_BUCK = STRESS_BUCK + local_STRESS_BUCK(i,:,:) 
    STRESS_EWALDD = STRESS_EWALDD + local_STRESS_EWALDD(i,:,:) 
    sum_buck = sum_buck + local_sum_buck(i) 
    sum_kin = sum_kin + local_sum_kin(i) 
    STRESS_KINETIC = STRESS_KINETIC + local_STRESS_KINETIC(i,:,:) 
end do 

скаляры и STRESS_KINETIC значения являются правильными, но STRESS_EWALDD неправильно, и я не могу понять, почему. Я пока не знаю о силах. Так что я очень благодарен за хит. Спасибо,

 Éric. 

Answer 1

Вы предприняли несколько неординарный подход к использованию OpenMP.

OpenMP обеспечивает положение reduction(op:vars), который выполняет редукцию op над локальными значениями переменных в vars и вы должны использовать его для STRESS_EWALD, sum_kin, sum_real и STRESS_KINETIC. Накопление силы для i-й атом должен работать, так как атомы в списке Verlet POINT заказываются (вы взяли код, который строит его из книги от Allen & Tildesley, справа?), Но не для j -й атом , Вот почему вы должны выполнять сокращение на них тоже. Не беспокойтесь, если вы прочитали в каком-то устаревшем руководстве OpenMP, что переменные сокращения должны быть скалярными - более новые версии OpenMP поддерживают сокращение по переменным массива в Fortran 90+.

!$OMP PARALLEL DO PRIVATE(....) 
!$OMP& REDUCTION(+:FDX,FDY,FDZ,sum_kin,sum_real,STRESS_EWALDD,STRESS_KINETIC) 
!$OMP& SCHEDULE(DYNAMIC) 
do i=1, atom_number  ! sum over atoms i 
type1 = ATOMS(i) 
do nj=POINT(i), POINT(i+1)-1 ! sum over atoms j of i's atoms list 
    j = LIST(nj) 
    type2 = ATOMS(j) 
    call get_ff_param(type1,type2,Aab,Bab,Cab,Dab) 
    call distance_sqr2(i,j,r,VECT_R) 
    call gettensor_rij(i,j,T) 
    r = sqrt(r) 
    if (r<coub_cutoff) then 
     sum_real = sum_real + (erfc(alpha_ewald*r) 
     derivative = -(erfc(alpha_ewald*r) 
     ff = - (1.0d0/r)*derivative 
     STRESS_EWALDD = STRESS_EWALDD + ff * T ! A 3x3 matrix  
     FDX(i) = FDX(i) + VECT_R(1) * ff 
     FDY(i) = FDY(i) + VECT_R(2) * ff 
     FDZ(i) = FDZ(i) + VECT_R(3) * ff 
     FDX(j) = FDX(j) - VECT_R(1) * ff 
     FDY(j) = FDY(j) - VECT_R(2) * ff 
     FDZ(j) = FDZ(j) - VECT_R(3) * ff  
    end if         
end do    ! sum over atoms j 

sum_kin = sum_kin + ATMMASS(i) * (VX(i)**2 + VY(i)**2 + VZ(i)**2) 
call gettensor_v(i,Q) 
STRESS_KINETIC = STRESS_KINETIC + ATMMASS(i) * Q 

end do 
!$OMP END PARALLEL DO 

Использование динамического цикла планирования уменьшит дисбаланс нагрузки, когда число соседей каждого атома будет дико.