13章

実際どのような分子とか結晶ができるのか、コンピュータで計算してみれば面白いけど難しいです。ポテンシャルに谷がたくさんあるので。
粒子間に働く力は距離の二乗に反比例しますが、ある粒子から距離Rだけ離れた位置にある粒子の数は距離の二乗に比例します。遠くに離れた粒子の影響をいつまでも受けてしまいます。たとえば力が遠くまで及ばない、分子同士の計算なら大きな立方体のかたまりを計算する時に表面がどうなっていても内部の計算には影響がありませんが、距離の二乗に反比例する力の場合は表面の状態が影響します。我々の世界でイオン結晶の計算をする時も同じ問題があります。

物理屋向けメモ：Yukawa型でmassに虚数をかける f(r)=exp(ar)/rのとき △ f = a^2 f

重力で引き合う3つの物体の軌跡

ネレオの力で引き合う3つの物体の軌跡

重力よりヘロヘロしてます

コード

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 3
#define DIM 3
double pos[N][DIM], pos_old[N][DIM], force[N][DIM];
double m=1.0;

void upd_verlet(double dt)
{
    int i,j;
    for(i=0;i<N;i++)
    {
	for(j=0;j<3;j++)
	{
	    double pos_new = 2*pos[i][j]-pos_old[i][j]+dt*dt/m*force[i][j];
	    pos_old[i][j] = pos[i][j];
	    pos[i][j] = pos_new;
	}
    }
}

// V= -cos(r)/r
// f= sin(r)/r +cos(r)/r^2

void cal_force1(double r[DIM], double f[DIM])
{
    int i;
    double r2=0;
    for(i=0;i<DIM;i++)r2 += r[i]*r[i];
    double r1=sqrt(r2);
#ifdef GRAVITY
    double f1= 1.0/r2;
#else
    double f1= ( sin(r1)/r1 +cos(r1)/r2);
#endif
    for(i=0;i<DIM;i++) f[i] = -f1 * r[i]/r1;
}

void cal_force()
{
    int i,j,k;
    for(i=0;i<N;i++)
    for(k=0;k<DIM;k++) force[i][k]=0.0;

    for(i=0;i<N;i++){
    for(j=0;j<i;j++){
	double rij[DIM], fij[DIM];
        for(k=0;k<DIM;k++) rij[k]= pos[i][k]-pos[j][k];
	cal_force1(rij, fij);

        for(k=0;k<DIM;k++) 
	{
	    force[i][k] += fij[k];
	    force[j][k] -= fij[k];
	}
    }}//ij
}

int main(int argc, char **argv)
{
    //random init
    int i,j,k;
    double v[N][DIM], vall[DIM];
    for(k=0;k<DIM;k++){
	vall[k]=0;
    for(i=0;i<N;i++){
	v[i][k] = 1.0*(rand()*1.0/RAND_MAX-0.5);
	vall[k] += v[i][k];
        }
        vall[k]/=N;
        for(i=0;i<N;i++) v[i][k] -= vall[k];
    }

    double dt=1.0e-4;
    for(i=0;i<N;i++){
    for(k=0;k<DIM;k++){
	pos[i][k]=pos_old[i][k]= rand()*30.0/RAND_MAX;
	pos_old[i][k] -= v[i][k]*dt;
    }}


    int t;
    for(t=0;t<3000000;t++)
    {
        cal_force();
        upd_verlet(dt);
	if(t%1000==0)
	{
            for(i=0;i<N;i++)
		printf("%f %f %f i%d\n", pos[i][0], pos[i][1], pos[i][2], i);

	}
    }
}