#include <stdio.h>
#include <math.h>
#include "mpi.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 N1 13               //各コアの担当する
#define N2 29               //最後の原子の番号
#define N3 49
#define N 100               //分子数
#define MD 2.29336758e24      //0.5/M_XENON

void input(FILE *file1,double r[][3],double v[][3]);
void next_r(double r[][3],double v[][3],double f[][3]);
double force(double f[][3],double r[][3],int my_rank);
void next_v(double v[][3],double total_f0[][3],double total_f1[][3],
	    double f1[][3],double total_potential,int n);
void output(FILE *file2,double r[][3]);

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

/***４コアで力を求める***/
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_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];
	}
      }
    }
  }
  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_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];
	}
      }
    }
  }
  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_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];
	}
      }
    }
  }
  else if(my_rank==3){
    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_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];
	}
      }
    }
  }
  return potential;
}

/***次ステップの位置を求める***/
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_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 total_f0[][3],double total_f1[][3],
	    double f1[][3],double total_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*total_f0[i][j]+DT*MD*total_f1[i][j];
      //力を初期化
      total_f0[i][j]=total_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;
  potential2*=4*E_XENON;
  //表示
  printf("%d:\t %le\t %le\t %le\n",n+1,kinetic,total_potential,kinetic+total_potential);
}

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