#include <stdio.h>
#include <math.h>
#define M_XENON 2.1802e-25  //キセノンの質量
#define E_XENON 3.20e-21    //キセノンのイプシロン
#define S_XENON 3.98e-10    //キセノンのシグマ
#define L_XENON 1.948e-9    //キセノンのセルの長さ
#define DT 1.0e-15          //微小時間
#define STEP 50000          //ステップ数
#define N 100               //分子数
#define MD 2.29336758e24     //0.5/M_XENON

void input(FILE *file1,double r[][3],double v[][3]);
void force(double f[][3],double r[][3]);
void next_r(double r[][3],double v[][3],double f0[][3]);
void next_v(double v[][3],double f0[][3],double f1[][3],int n);
void output(FILE *file2,double r[][3]);
double potential,kinetic;

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

  file1=fopen("xenon.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<STEP;n++){  
    /***次ステップの位置を求める***/
    next_r(r,v,f0);
    /***次ステップの力を求める***/
    force(f1,r);
    /***次ステップの速度を求め、エネルギーの表示***/
    next_v(v,f0,f1,n);
    /***50ステップごとに座標を書き込み***/
    if(n%50==0)
       output(file2,r);
  }
  fclose(file1);
  fclose(file2);
  return 0;
}
/**************************************************************************************************************/
/***初期位置と初期速度を取り込み***/
void input(FILE *file1,double r[][3],double v[][3]){
  int i,j;
  //位置
  for(i=0;i<N;i++)
    for(j=0;j<3;j++)
      fscanf(file1,"%le",&r[i][j]);
  //速度
  for(i=0;i<N;i++)
    for(j=0;j<3;j++)
      fscanf(file1,"%lf",&v[i][j]);
}

/***力を求める***/
void force(double f[][3],double r[][3]){
  int i,j,k;
  double dis;  //２分子間の距離
  double vec[3];  //分子iから分子jへのベクトル
  double t1,t2,t3,t4,t5;
  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_XENON<vec[k])
	  vec[k]-=L_XENON;
	else if(-0.5*L_XENON>vec[k])
	  vec[k]+=L_XENON;
      }
      dis=pow(vec[0],2)+pow(vec[1],2)+pow(vec[2],2);
      t1=pow(S_XENON,2)/dis;
      t2=t1*t1*t1;
      t3=pow(t2,2);
      t4=1/dis;
      potential+=t3-t2;
      t5=24*E_XENON*t4*(2*t3-t2);
      //力を求める（作用反作用の法則）
      for(k=0;k<3;k++){
	f[i][k]+=t5*vec[k];
	f[j][k]-=t5*vec[k];
      }
    }
  }
}

/***次ステップの位置を求める***/
void next_r(double r[][3],double v[][3],double f0[][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*f0[i][j];
      //周期境界条件
      if(r[i][j]>0.5*L_XENON)
	r[i][j]-=L_XENON;
      if(r[i][j]<-0.5*L_XENON)
	r[i][j]+=L_XENON;
    }
}

/***次ステップの速度を求め、エネルギーの表示***/
void next_v(double v[][3],double f0[][3],double f1[][3],int n){
  int i,j;
  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][0],2)+pow(v[i][1],2)+pow(v[i][2],2);
  }
  kinetic*=0.5*M_XENON;
  potential*=4*E_XENON;
  //表示
  printf("%d:\t %le\t %le\t %le\n",n+1,kinetic,potential,kinetic+potential);
}

/***座標を書き込み***/
void output(FILE *file2,double r[][3]){
  int i;
  for(i=0;i<N;i++)
    //位置を10^10倍する
    fprintf(file2,"%d Xe %f %f %f\n",i+1,r[i][0]*100000*100000,r[i][1]*100000*100000,r[i][2]*100000*100000);
}