如何用OpenMP优化小系统经典分子动力学并行模拟

我尝试用OpenMP并行化经典分子动力学程序，但核心困境是需要串行运行多个**小系统（原子数<100）**的模拟，而非单个大系统。以13个原子的系统为例，运行50次各含约1000步的模拟后，并行化反而导致性能大幅下降：

• 无OpenMP：耗时0.15秒
• 使用omp parallel for reduction(+:E)：耗时0.50秒
• 使用omp parallel for shared(E)：耗时1.50秒
  并行化后最慢速度甚至比串行慢了10倍！

原始程序结构 #

最初我对原子循环进行并行化，核心函数如下：

double evalLJQs_ng4_omp( const Quat4i* neighs, double Rdamp=1.0 ){
    double R2damp = Rdamp*Rdamp;
    double E =0;
    #pragma omp parallel for reduction(+:E)
    for(int i=0; i<natoms; i++){
       E += evalLJQs_ng4_atom( i, neighs, R2damp );
    }
    return E;
}

循环中调用evalLJQs_ng4_atom函数，该函数读取全局数组，计算Lennard-Jones和库仑相互作用能与力，最终将力写入全局数组fs的对应位置（无写入冲突）：

double evalLJQs_ng4_atom( int ia, const Quat4i* neighs, double R2damp ){
    Vec3d        fi   = Vec3dZero;
    Vec3d        pi   = ps    [ia];           // 全局读取 ps[ia]
    const Vec3d  REQi = REQs  [ia];           // 全局读取 REQs[ia]
    const Quat4i ngs  = neighs[ia];           // 全局读取 neighs[ia]
    double         E  = 0;
    for(int j=0; j<natoms; j++){
        if(ia==j) continue;                                             // 跳过自相互作用
        if( (ngs.x==j)||(ngs.y==j)||(ngs.z==j)||(ngs.w==j) ) continue;  // 跳过成键邻居
        Vec3d fij   = Vec3dZero;
        const Vec3d pj    = ps  [j];                        // 全局读取 ps[j]
        const Vec3d REQj  = REQs[j];                        // 全局读取 REQs[j]
        Vec3d REQij; combineREQ( REQj, REQi, REQij );       // 纯函数（无全局变量），L-J混合规则
        E          += getLJQ( pj-pi, REQij, R2damp, fij );  // 纯函数（无全局变量），计算L-J与库仑相互作用
        fi.add(fij);
    }
    fs[ia].add(fi);           // 全局写入 fs[ia]，无冲突
    return E;
}

关键说明 #

1. 上述#pragma omp parallel是程序中唯一的OpenMP并行区域
2. 该并行函数仅为整个分子动力学循环的一小部分，其他模块（如Verlet传播器）未并行——小系统下这些模块并行收益极低
3. 推测性能下降的核心原因是每步模拟重复初始化/销毁并行上下文的开销，希望仅通过#pragma指令进行非侵入式优化，避免大幅修改程序结构

优化尝试 #

• 替换shared(E)为#pragma omp parallel for reduction(+:E)，减少全局变量的竞争开销
• 将并行区域移至整个动力学循环之外，避免每步重复创建线程上下文，最终得到可用的优化版本：

int run_omp( int niter, double dt, double Fconv, double Flim ){
    double F2conv = Fconv*Fconv;
    double E=0,F2=0;
    int    itr=0;
    #pragma omp parallel shared(E,F2) private(itr)
    for(itr=0; itr<niter; itr++){
        // 仅由单个线程执行初始化
        #pragma omp single
        {E=0;F2=0;}
        // ------ 计算MMFF能量
        #pragma omp for reduction(+:E)
        for(int ia=0; ia<natoms; ia++){ 
            if(verbosity>3)printf( "atom[%i]@cpu[%i/%i]\n", ia, omp_get_thread_num(), omp_get_num_threads() );
            if(ia<nnode)E += eval_atom(ia);
            E += evalLJQs_ng4_PBC_atom( ia ); 
        }
        // ---- 组装结果（需等待所有原子计算完成）
        #pragma omp for
        for(int ia=0; ia<natoms; ia++){
            assemble_atom( ia );
        }
        // ------ 原子运动更新
        #pragma omp for reduction(+:F2)
        for(int i=0; i<nvecs; i++){
            F2 += move_atom_MD( i, dt, Flim, 0.99 );
        }
        #pragma omp single 
        if(verbosity>2){printf( "step[%i] E %g |F| %g ncpu[%i] \n", itr, E, F2, omp_get_num_threads() );}
    }
    return itr;
}

内容的提问来源于stack exchange，提问作者Prokop Hapala