#include <stdio.h>
#include <math.h>
#include "mpi.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 N1 13               //各コアの担当する
#define N2 30               //最後の原子の番号
#define N3 50
#define N 100               //分子数
#define MD 7.537499e24      //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 total_f0[][3],double total_f1[][3],
	    double f1[][3],int n,int my_rank);
double force(double f[][3],double r[][3],int my_rank);

int main(int argc,char *argv[]){
  FILE *file1,*file2;
  double r[N][3];   	   //位置
  double v[N][3];   	   //速度
  double f0[N][3]={{0}};   //今の力
  double total_f0[N][3];   //f0の和
  double f1[N+1][3]={{0}}; //次の力
  double total_f1[N+1][3]; //f1の和
  int my_rank,nproc;
  int n;
  
  /*MPI初期化*/
  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
  MPI_Comm_size(MPI_COMM_WORLD,&nproc);
  /*初期処理*/
  file1=fopen("argon.inp","r");
  input(file1,r,v);
  if(my_rank==1){
    printf("cycle \t kinetic \t potential \t total\n");
    file2=fopen("mol.arc","w");
  }
  /*初期の力とその和*/
  force(f0,r,my_rank);
  MPI_Allreduce(f0,total_f0,3*N,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
  /*以下を繰り返す*/
  for(n=0;n<STEP;n++){  
    /*次の位置、力、ポテンシャルとその和*/
    next_r(r,v,total_f0);
    //f1[N][0]をポテンシャルの入れ物にする
    f1[N][0]=force(f1,r,my_rank);
    MPI_Allreduce(f1,total_f1,3*N+1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    /*次の速度、力の初期化、エネルギーの計算と表示*/
    result(v,total_f0,total_f1,f1,n,my_rank);
    /*50ステップごとに座標を書き込み*/
    if(n%50==0&&my_rank==1)
       output(file2,r);
  }
  /*終了処理*/
  fclose(file1);
  if(my_rank==1)
    fclose(file2);
  MPI_Finalize();
  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;
      else if(r[i][j]<-0.5*L_ARGON)
	r[i][j]+=L_ARGON;
    }
}

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


/*力*/
double force(double f[][3],double r[][3],int my_rank){
  int i,j,k;
  double dis;     //距離の２乗
  double vec[3];  //分子iから分子jへのベクトル  
  double t1,t2,t3,t4,t5;
  double potential=0;
  if(my_rank==0){
    for(i=0;i<N1;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];
	}
      }
    }
  }
  else if(my_rank==1){
    for(i=N1;i<N2;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];
	}
      }
    }
  }
  else if(my_rank==2){
    for(i=N2;i<N3;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];
	}
      }
    }
  }
  else{
    for(i=N3;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;
}