5

Pythonでアルゴン液体の分子動力学シミュレーションを行っています。安定バージョンを実行していますが、100アトムを超えると実行速度が遅くなります。ボトルネックは、次のネストされたforループであると特定しました。これは、main.pyスクリプトから呼び出される関数内に配置された力の計算です。

def computeForce(currentPositions):
    potentialEnergy = 0
    force = zeros((NUMBER_PARTICLES,3))
    for iParticle in range(0,NUMBER_PARTICLES-1):
        for jParticle in range(iParticle + 1, NUMBER_PARTICLES):
            distance = currentPositions[iParticle] - currentPositions[jParticle]
            distance = distance - BOX_LENGTH * (distance/BOX_LENGTH).round()
            #note: this is so much faster than scipy.dot()
            distanceSquared = distance[0]*distance[0] + distance[1]*distance[1] + distance[2]*distance[2]            
            if distanceSquared < CUT_OFF_RADIUS_SQUARED:
                r2i = 1. / distanceSquared
                r6i = r2i*r2i*r2i
                lennardJones = 48. * r2i * r6i * (r6i - 0.5)
                force[iParticle] += lennardJones*distance
                force[jParticle] -= lennardJones*distance
                potentialEnergy += 4.* r6i * (r6i - 1.) - CUT_OFF_ENERGY
return(force,potentialEnergy)

大文字の変数は定数であり、config.pyファイルで定義されています。「currentPositions」は、粒子数が3の行列です。

ネストされたforループのバージョンをscipy.weaveで既に実装しました。これは、次のWebサイトから着想を得ています：http ：//www.scipy.org/PerformancePython 。

しかし、私は柔軟性の喪失が好きではありません。これをforループで「ベクトル化」することに興味があります。私はそれがどのように機能するのか本当にわかりません。誰かが私にこれを教える手がかりや良いチュートリアルを教えてもらえますか？

4

3 に答える 3

4

以下は私のコードのベクトル化されたバージョンです。1000ポイントのデータセットの場合、私のコードは元のコードよりも約50倍高速です。

In [89]: xyz = 30 * np.random.uniform(size=(1000, 3))

In [90]: %timeit a0, b0 = computeForce(xyz)
1 loops, best of 3: 7.61 s per loop

In [91]: %timeit a, b = computeForceVector(xyz)
10 loops, best of 3: 139 ms per loop

コード：

from numpy import zeros

NUMBER_PARTICLES = 1000
BOX_LENGTH = 100
CUT_OFF_ENERGY = 1
CUT_OFF_RADIUS_SQUARED = 100

def computeForceVector(currentPositions):
    potentialEnergy = 0
    force = zeros((NUMBER_PARTICLES, 3))
    for iParticle in range(0, NUMBER_PARTICLES - 1):
        positionsJ =  currentPositions[iParticle + 1:, :]
        distance = currentPositions[iParticle, :] - positionsJ
        distance = distance - BOX_LENGTH * (distance / BOX_LENGTH).round()
        distanceSquared = (distance**2).sum(axis=1)
        ind = distanceSquared < CUT_OFF_RADIUS_SQUARED

        if ind.any():
            r2i = 1. / distanceSquared[ind]
            r6i = r2i * r2i * r2i
            lennardJones = 48. * r2i * r6i * (r6i - 0.5)
            ljdist = lennardJones[:, None] * distance[ind, :]
            force[iParticle, :] += (ljdist).sum(axis=0)
            force[iParticle+1:, :][ind, :] -= ljdist
            potentialEnergy += (4.* r6i * (r6i - 1.) - CUT_OFF_ENERGY).sum()
    return (force, potentialEnergy)

コードが同じ結果を生成することも確認しました

于 2013-02-13T15:35:40.590 に答える
3

純粋なPythonでMDエンジンのようなものを書くのは遅くなります。Numba（ http://numba.pydata.org/）またはCython（http://cython.org/ ）のいずれかを見てみます。Cython側では、cythonを使用した単純なランジュバン動力学エンジンを作成しました。これは、開始するための例として役立つ可能性があります。

https://bitbucket.org/joshua.adelman/pylangevin-integrator

私がとても気に入っているもう1つのオプションは、OpenMMを使用することです。MDエンジンのすべての部分をまとめたり、カスタムフォースを実装したりできる、Pythonラッパーがあります。GPUデバイスをターゲットにする機能もあります。

https://simtk.org/home/openmm

ただし、一般的に、高度に調整されたMDコードは非常に多くあるため、ある種の一般的な教育目的でこれを実行しない限り、独自のコードを最初から作成することは意味がありません。いくつか例を挙げると、主要なコードのいくつか：

于 2013-02-13T14:34:03.393 に答える
1

この投稿を完成させるために、Cコードで織り込まれた実装を追加します。これを機能させるには、ウィーブとコンバーターをインポートする必要があることに注意してください。さらに、weaveは現在python2.7でのみ機能します。すべての助けにもう一度感謝します！これは、ベクトル化されたバージョンよりもさらに10倍高速に実行されます。

from scipy import weave
from scipy.weave import converters
def computeForceC(currentPositions):        
    code = """
    using namespace blitz;
    Array<double,1> distance(3);
    double distanceSquared, r2i, r6i, lennardJones;
    double potentialEnergy = 0.;

    for( int iParticle = 0; iParticle < (NUMBER_PARTICLES - 1); iParticle++){
        for( int jParticle = iParticle + 1; jParticle < NUMBER_PARTICLES; jParticle++){
            distance(0) = currentPositions(iParticle,0)-currentPositions(jParticle,0);
            distance(0) = distance(0) - BOX_LENGTH * round(distance(0)/BOX_LENGTH);
            distance(1) = currentPositions(iParticle,1)-currentPositions(jParticle,1);
            distance(1) = distance(1) - BOX_LENGTH * round(distance(1)/BOX_LENGTH);
            distance(2) = currentPositions(iParticle,2)-currentPositions(jParticle,2);
            distance(2) = distance(2) - BOX_LENGTH * round(distance(2)/BOX_LENGTH);
            distanceSquared = distance(0)*distance(0) + distance(1)*distance(1) + distance(2)*distance(2);
            if(distanceSquared < CUT_OFF_RADIUS_SQUARED){
                r2i = 1./distanceSquared;
                r6i = r2i * r2i * r2i;
                lennardJones = 48. * r2i * r6i * (r6i - 0.5);
                force(iParticle,0) += lennardJones*distance(0);
                force(iParticle,1) += lennardJones*distance(1);
                force(iParticle,2) += lennardJones*distance(2);
                force(jParticle,0) -= lennardJones*distance(0);
                force(jParticle,1) -= lennardJones*distance(1);
                force(jParticle,2) -= lennardJones*distance(2);
                potentialEnergy += 4.* r6i * (r6i - 1.)-CUT_OFF_ENERGY;

                }

            }//end inner for loop
    }//end outer for loop
    return_val = potentialEnergy;

    """
    #args that are passed into weave.inline and created inside computeForce
    #potentialEnergy = 0.
    force = zeros((NUMBER_PARTICLES,3))

    #all args
    arguments = ['currentPositions','force','NUMBER_PARTICLES','CUT_OFF_RADIUS_SQUARED','BOX_LENGTH','CUT_OFF_ENERGY']
    #evaluate stuff in code
    potentialEnergy = weave.inline(code,arguments,type_converters = converters.blitz,compiler = 'gcc')    

    return force, potentialEnergy
于 2013-02-15T08:13:08.833 に答える