使用 Fortran 90 的 OpenMP 和数组求和

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 do”指令并行化这个双循环，但遇到了张量 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

当每个原子的邻居数量变化很大时，使用动态循环调度将减少负载不平衡。