0

我正在尝试用 C/C++ 编写一个简单的 MD 程序（我习惯于 C 但我正在尝试学习 C++，所以我的代码有点“混合”......我知道这是次优的，我一旦我完全理解它就会转向完整的C++）。

一切似乎都在运行，但我在动能方面存在分歧，系统不会热化，温度（prop to K）一步即可从阶（10°K）变为阶（10000°K）。

我正在使用 0.002 的低时间步长（模拟总时间：30），所以我不应该有这个巨大的错误......

这是我的代码，如果不清楚，我可以尝试更好地解释它

int main(){

...

int n, t, m, i;
double r, K, U, E,P,  totalE, temperature, d, x,y,z, temp;

...

double data[5][PARTICELLE], vel[3][PARTICELLE],dataNew[3][PARTICELLE]; //0,1,2 are x,y,z. 3, 4 for data are Energy and Pressure
double force[3][PARTICELLE], forceNew[3][PARTICELLE];
double velQ[PARTICELLE]; //square velocity

ofstream out(OUTDATA);

//inizio MD

for(t=0; t<PASSI; t++){
    //inizialization
    K=0;
    U=0;
    E=0;
    P=0;
    fill(data[3], data[3]+PARTICELLE, 0); //E=0 for each particle
    fill(data[4], data[4]+PARTICELLE, 0);
    fill(velQ, velQ+PARTICELLE, 0);
    for(i=0; i<3; i++){
        fill(force[i], force[i]+PARTICELLE, 0);
        fill(forceNew[i], forceNew[i]+PARTICELLE, 0);
    }
    for(n=0; n<PARTICELLE; n++) { //for on the n_ particle. A step is a move of n=PARTICELLE particles
        for (i = 0; i < 3; i++) { //compute vSquare
            velQ[n] += vel[i][n] * vel[i][n];
        }
        K += 0.5 * MASSA * velQ[n]; //compute Kinetic Energy
        for(m=0; m<PARTICELLE; m++){ //loop on m!=n to compute F, E, P
            if(m!=n){
                r=0;
                for(i=0; i<3; i++){ //calculation of radius and x,y,z 
                    d = data[i][m] - data[i][n];
                    d = d - (NINT(d / LATO) * LATO);
                    if(i==0)x=d;
                    if(i==1)y=d;
                    if(i==2)z=d;
                    r += d * d;
                }
                //if (t<2) cout << "x y z" << x << " " << y << " " << z << endl;
                r=sqrt(r);
                if (r < R) {
                    data[3][n] += energy(r); //update Energy of n
                    for(i=0; i<3; i++){
                        if(i==0)temp=x;
                        if(i==1)temp=y;
                        if(i==2)temp=z;
                        force[i][n]+=forza(r,temp); //compute force (cartesian components)
                        //if(t<2)cout << "force " <<n << " " << m << " "<< force[i][n] << endl;
                    }
                    if (m < n)data[4][n] += (-energy(r) * (1 + r)); //pressure
                }
            }
        }
        U+=data[3][n]; //total potential energy 
        P+=data[4][n]; //total pressure
        for (i = 0; i < 3; i++) { //Verlet update, part 1
            dataNew[i][n] = data[i][n] + vel[i][n] * DeltaT + 0.5 * force[i][n] * DeltaT * DeltaT / MASSA;
        }
        for(m=0; m<PARTICELLE; m++){ //update force
            if(m!=n){
                r=0;
                for(i=0; i<3; i++){
                    d = data[i][m] - dataNew[i][n];
                    d = d - (NINT(d / LATO) * LATO);
                    if(i==0)x=d;
                    if(i==1)y=d;
                    if(i==2)z=d;
                    r += d * d;
                }
                r=sqrt(r);
                if (r < R) {
                    for(i=0; i<3; i++) {
                        if (i == 0)temp = x;
                        if (i == 1)temp = y;
                        if (i == 2)temp = z;
                        forceNew[i][n] += forza(r, temp);

                    }
                }
            }
        }
        for(i=0; i<3; i++){ //new position and Verlet part 2 
            data[i][n]=dataNew[i][n];
            vel[i][n]=vel[i][n] + DeltaT * 0.5*(forceNew[i][n] + force[i][n]) / MASSA;
        }
    }
    totalE=U+K; //total energy
    temperature = 2*K/(PARTICELLE*3);
    out << t*DeltaT << " " << U  << " " << P  << " " << totalE  << " " << temperature << endl;
}
out.close();
return 0;
}

我的系统在潜在的 e^-r/r 下，所以我有：

double energy( double r){
return (A*SIGMA*exp(-r/SIGMA)/r);
}
double forza(double r, double h){ //h is for x,y,z
    double bubba;
    bubba= (A*SIGMA*(exp(-r)*h*(r+1)/(r*r*r)));
    return bubba;
}

谢谢你的帮助。自四月以来我一直在研究此代码，但我仍然没有解决方案......

编辑：更清楚：大写术语 和DeltaT是定义的值DEFINE

4

0 回答 0