如何循环潜力？答案

【问题标题】：How to loop the potential?如何循环潜力？
【发布时间】：2015-10-26 23:07:16
【问题描述】：

我目前正在研究溶液中聚合物的分子动力学模拟，这是计算系统势能和施加在每个单体上的力的子程序之一。

function eval_force()
% eval_force.m IS USED FOR EVALUATING FORCE


% THE STRATEGY USUALLY ADOPTED FOR A LENNARD-JONES OR AS A MATTER OF FACT
% ANY PAIR-WISE INTERACTING SYSTEM IS AS FOLLOWS:
% 1. EVALUATE THE DISTANCE BETWEEN TWO PAIRS OF ATOMS
% 2. ENSURE THAT MINIMUM IMAGE CONVENTION (MIC) IS FOLLOWED
% 3. IF THE DISTANCE OBTAINED THROUGH MIC IS GREATER THAN THE CUT OFF
% DISTANCE MOVE TO NEXT PAIR
% 4. ELSE EVALUATE POTENTIAL ENERGY AND CALCULATE FORCE COMPONENTS
% 5. F(i,j) = -F(j,i)

global MASS KB TEMPERATURE NUM_ATOMS LENGTH TSTEP;
global EPS SIG R_CUT GAMMA POT_E;
global POSITION VELOCITY FORCE STO;

dr = zeros(3,1);
drh = zeros(3,1);
FORCE(:) = 0.0;
POT_E = 0.0;

for ( i=1:NUM_ATOMS )
    for ( j=i+1:NUM_ATOMS )

        dist2 = 0.0; % VARIABLE dist2 STORES DISTANCE BETWEEN PAIR (i,j)

        % FIRST FIND OUT THE DIFFERENCE IN X,Y AND Z COORDINATES
        % VARIABLE dr IS USED FOR THIS PURPOSE
        for(k = 1:3)
            dr(k) = POSITION(i,k) - POSITION(j,k);

            % THESE STEPS ENSURE MINIMUM IMAGE CONVENTION IS FOLLOWED
            if(dr(k) > LENGTH/2.0)
                dr(k) = dr(k) - LENGTH;
            end
            if(dr(k) < -LENGTH/2.0)
                dr(k) = dr(k) + LENGTH;
            end
            % MINIMUM IMAGE CONVENTION ENDS HERE

            dist2 = dist2 + dr(k)*dr(k); % dist2 IS BASED UPON MIC
        end

        if(dist2 <= R_CUT*R_CUT) % IF THE CUT OFF CRITERIA IS SATISFIED
            dist2i = power(SIG,2)/dist2;
            dist6i = power(dist2i,3);
            dist12i = power(dist6i,2);
            POT_E = POT_E + EPS * (dist12i - 2*dist6i) + 33.34 * EPS *     power(sqrt(dist2) - SIG,2)/(2 * power(SIG,2)); % STORES THE POTENTIAL ENERGY

            Ff = 12.0 * EPS * (dist12i-dist6i) - 33.34 * EPS * (sqrt(dist2) - SIG)/(dist2i * sqrt(dist2) * power(SIG,2));
            Ff = Ff * dist2i;

            for(k = 1:3)
                FORCE(i,k) = FORCE(i,k) + Ff*dr(k)- GAMMA*VELOCITY(i,k);
                FORCE(j,k) = FORCE(j,k) - Ff*dr(k)- GAMMA*VELOCITY(j,k);
            end
        end
   end
end

end

如何为 POT_E 下的“33.34 * EPS * power(sqrt(dist2) - SIG,2)/(2 * power(SIG,2))”部分创建一个循环，这是谐波电位，所以它只计算最近原子之间的距离（对于 j=i+1 到 4）。

【问题讨论】：

您知道，与 matlab 相比，在低级语言中（例如 fortran，该代码最初的意图是，至少表面上是这样）您将获得巨大的加速？对于任何严肃的应用程序，我建议在完成原型设计和概念验证后迁移到较低级别的语言。

标签： matlab simulation

【解决方案1】：

在快速浏览您的代码后，我建议以下想法：

使用pdist 一次计算所有距离并将其存储在矩阵中，例如AllDist。
基于LENGTH的条件可以直接应用于AllDist。
由于您需要找到四个最近的邻居，您需要有一个循环遍历AllDist 的行（或列），它对当前行（或列）进行排序并为您提供四个最近的邻居。请注意，对于每个原子，您将获得 0 作为最近距离，因为这是“自距离”。忽略这个。
如果您可以访问 Matlab 的并行计算工具箱，请尝试在适当的地方使用它 (parfor) 来加速您的模拟。

【讨论】：