#include <stdio.h>
#include <math.h>
#define M_ARGON 6.6335e-26  //アルゴンの質量
#define E_ARGON 1.67e-21    //アルゴンのイプシロン
#define S_ARGON 3.40e-10    //アルゴンのシグマ
#define L_ARGON 1.680e-9    //アルゴンのセルの長さ
#define DT 1.0e-15          //微小時間
#define STEP 50000          //ステップ数
#define N 100               //分子数
#define MD 7.5374990578e24  //0.5/M_ARGON

void input(FILE *file,double r[][3],double v[][3]);
void output(FILE *file,double r[][3]);
void next_r(double r[][3],double v[][3],double f[][3]);
void result(double v[][3],double f0[][3],double f1[][3],double potential,int n);
double force(double f[][3],double r[][3]);

int main(void){
  FILE *file1,*file2;
  double r[N][3];   	//位置
  double v[N][3];   	//速度
  double f0[N][3]={{0}};//今の力
  double f1[N][3]={{0}};//次の力
  double potential;
  int n;

  /*初期処理*/
  file1=fopen("argon.inp","r");
  file2=fopen("mol.arc","w");
  fprintf(file2,"100\n");
  input(file1,r,v);
  printf("cycle \t kinetic \t potential \t total\n");
  /*以下を繰り返す*/
  for(force(f0,r),n=0;n<STEP;n++){  
    /*次の位置、力、ポテンシャル*/
    next_r(r,v,f0);
    potential=force(f1,r);
    /*次の速度、力の初期化、エネルギーの計算と表示*/
    result(v,f0,f1,potential,n);
    /*50ステップごとに座標を書き込み*/
    if(n%50==0)
       output(file2,r);
  }
  /*終了処理*/
  fclose(file1);
  fclose(file2);
  return 0;
}

/*初期位置と初期速度を取り込み*/
void input(FILE *file,double r[][3],double v[][3]){
  int i,j;
  //位置
  for(i=0;i<N;i++)
    for(j=0;j<3;j++)
      fscanf(file,"%le",&r[i][j]);
  //速度
  for(i=0;i<N;i++)
    for(j=0;j<3;j++)
      fscanf(file,"%lf",&v[i][j]);
}

/*座標を書き込み*/
void output(FILE *file,double r[][3]){
  int i;
  fprintf(file,"100\n");
  for(i=0;i<N;i++)
    fprintf(file,"%d Ar %f %f %f\n",i+1,
	    r[i][0]*100000*100000,r[i][1]*100000*100000,r[i][2]*100000*100000);
}

/*次の位置*/
void next_r(double r[][3],double v[][3],double f[][3]){ 
  int i,j;
  for(i=0;i<N;i++)
    for(j=0;j<3;j++){
      r[i][j]+=DT*v[i][j]+DT*DT*MD*f[i][j];
      //周期境界条件
      if(r[i][j]>0.5*L_ARGON)
	r[i][j]-=L_ARGON;
      if(r[i][j]<-0.5*L_ARGON)
	r[i][j]+=L_ARGON;
    }
}

/*次の速度、力の初期化、エネルギーの計算と表示*/
void result(double v[][3],double f0[][3],double f1[][3],double potential,int n){
  int i,j;
  double kinetic=0;
  for(i=0;i<N;i++){
    for(j=0;j<3;j++){
      //次の速度
      v[i][j]+=DT*MD*f0[i][j]+DT*MD*f1[i][j];
      //力を初期化
      f0[i][j]=f1[i][j];
      f1[i][j]=0;
      //エネルギーの計算と表示
      kinetic+=pow(v[i][j],2);
    }
  }
  kinetic*=0.5*M_ARGON;
  potential*=4*E_ARGON;
  printf("%d:\t %le\t %le\t %le\n",n+1,kinetic,potential,kinetic+potential);
}

/*力*/
double force(double f[][3],double r[][3]){
  int i,j,k;
  double dis;     //２分子間の距離
  double vec[3];  //分子iから分子jへのベクトル
  double t1,t2,t3,t4,t5;
  double potential=0;
  for(i=0;i<N;i++){
    for(j=i+1;j<N;j++){
      for(k=0;k<3;k++){
	vec[k]=r[i][k]-r[j][k];
	//最も近い分子を選ぶ
	if(0.5*L_ARGON<vec[k])
	  vec[k]-=L_ARGON;
	else if(-0.5*L_ARGON>vec[k])
	  vec[k]+=L_ARGON;
      }
      dis=pow(vec[0],2)+pow(vec[1],2)+pow(vec[2],2);
      t1=pow(S_ARGON,2)/dis;
      t2=t1*t1*t1;
      t3=1/dis;
      t4=pow(t2,2);
      //ポテンシャルを計算
      potential+=t4-t2;
      t5=24*E_ARGON*t3*(2*t4-t2);
      //力を計算
      for(k=0;k<3;k++){
	f[i][k]+=t5*vec[k];
	f[j][k]-=t5*vec[k];
      }
    }
  }
  return potential;
}