热力学分布用matlab,关于径向分布函数RDFs的代码(求指正)！

我把写的C代码(源自于《the art of MD)附上，自己参考一下，

void EvalRdf ()

{

VecR dr;

real deltaR, normFac, rr;

int j1, j2, n;

if (countRdf == 0)

{

for (n = 0; n < sizeHistRdf; n ++)

{

histRdf[n] = 0.;

}

}

deltaR = rangeRdf / sizeHistRdf;

for (j1 = 0; j1 < nMol - 1; j1 ++)

{

for (j2 = j1 + 1; j2 < nMol; j2 ++)

{

VSub (dr, mol[j1].r, mol[j2].r);

VWrapAll (dr);

rr = VLenSq (dr);

if (rr < Sqr (rangeRdf))

{

n = sqrt (rr) / deltaR;

++ histRdf[n];

}

}

}

++ countRdf;

if (countRdf == limitRdf)

{

normFac = VProd (region) / (2. * M_PI * Cube (deltaR) * Sqr (nMol) * countRdf);

for (n = 0; n < sizeHistRdf; n ++)

{

histRdf[n] *= normFac / Sqr (n - 0.5);

}

PrintRdf (stdout);

countRdf = 0;

}

}

fortan的一些结构写法我不太明白。你自己对照改一下就行了，算法都一样。