#include <stdio.h>
#include <fstream.h>

double rint(double x) {
  	return ( (x<0.)? -floor(-x+.5):floor(x+.5) );
}    //rint：四捨五入を与える関数の定義
		

int main(){
	
//=======(�@ 変数と出力ファイルの設定）================================

	const double dt      = 4.0;          // 積分時間 [fs]
	const int    maxstep = 10000;         // ステップ数
	const double Tset    = 700.;         // 温度 [K]
	const double Pset    = 300.*101.3e+3;// 圧力 [Pa]
	const int    N       = 50;          // 原子数 [個/u.c.]
	const int    spnum   = 1;            // 成分数
	
	const double KB      = 1.38065e-23;// ボルツマン定数[J/K]
	
	//LJポテンシャル変数 単位はsig：Å,eps：J
	//原子の種類ごとにパラメータを設定			
	double sig[spnum],eps[spnum];
	sig[0] = 3.4; 
	eps[0] = 121.0*KB;	 
	
	double cutoff  = 9.0;             //カットオフ [A]

	double w[N];                      // 原子の重さ [g/mol]
	double rx[N],ry[N],rz[N];         // 原子の座標 [A]
	double ax[N],ay[N],az[N];         // MSD用の変数
	double vx[N],vy[N],vz[N];         // 原子の速度 [m/s]
	double fx[N],fy[N],fz[N];         // 各原子の力
	int    sp[N];
	double cx,cy,cz;      // 単位セルの大きさ
	double T;             // 系の温度 [K]
	double KE,PE;         // エネルギーの保存用
	int    i,j,k,m;       // 整数型の変数
	double r,s,t,F,x,y,z; // 実数型の変数
	
    //出力用ファイルを作製する                               
	ofstream fpos("position.dat");// 座標出力用
	ofstream fene("energy.dat");  // エネルギー値出力用
	ofstream fmsd("msd.dat");     // 座標出力用

//=======(�A ＭＤ「単位セル」の大きさの計算）=========================

	double v;                       // 単位セルの体積 [A^3]
	const double AVO = 6.02e+23;    // アボガドロ数
	const double R   = 8.314;       // 気体定数 [J/(K mol)]
	
	v = ((double)N/AVO*R*Tset)/Pset;    // PV=nRTから体積を計算
	cx = cy = cz = pow(v,1./3)*1.0e+10; // 立方体として各辺を計算
	
	cout <<"セルの一辺の長さは、"<<cx<<" です。"<<endl;           

//=======(�B 原子の初期配置の計算）===================================

	double rmin = 2.;          //最小原子間距離 [A]
	
	for(i=0;i<N;i++){
		rx[i] = (double)rand()/(RAND_MAX)*cx; // ランダムに座標を発生
		ry[i] = (double)rand()/(RAND_MAX)*cy; // 単位は[Å]
		rz[i] = (double)rand()/(RAND_MAX)*cz;
		
		w[i]  = 39.948;           // 原子量の入力 [g/mol]
		sp[i] = 0;                // 原子種番号の入力

		ax[i] = 0.0;
		ay[i] = 0.0;
		az[i] = 0.0;
										
		for(j=0;j<i;j++){
			r  = (rx[j]-rx[i])*(rx[j]-rx[i])
			      + (ry[j]-ry[i])*(ry[j]-ry[i])
			       + (rz[j]-rz[i])*(rz[j]-rz[i]);
			if(r<rmin*rmin) {
			  cout<<i<<"番目の原子位置が重なったので再計算します。"
				        <<endl; i--;}
		}
	}
	cout <<N<<"個の原子のランダム配置が完了しました。"<<endl;	
	
//=======(�C 原子初期速度の設定）=====================================

	for(i=0;i<N;i++){
		vx[i] = 0.0;
		vy[i] = 0.0;
		vz[i] = 0.0;
	}
	cout<<"原子の初期速度を設定しました。"<<endl;	

//=======(�D 反復計算の開始）=========================================
		
	double dx,dy,dz,sigIJ,epsIJ,r2;  //実数型変数の定義
	
	for(k=0;k<maxstep;k++){    //ステップ(k)のループの始まり

		
//=======(�E 原子間力の計算）=========================================

		for(i=0;i<N;i++){           //力の配列の初期化
			fx[i]=fy[i]=fz[i]=0.0;
		}
		PE = 0.0;       //ポテンシャルエネルギー変数の初期化

		for(i=0;i<N;i++){           // 原子iについてのfor文
			for(j=0;j<i;j++){       // 原子ｊについてのfor文
				
				dx = rx[j]-rx[i];   // 原子間のｘベクトル計算
				dy = ry[j]-ry[i];   // 原子間のｙベクトル計算
				dz = rz[j]-rz[i];   // 原子間のｚベクトル計算
						
			  	dx = dx-rint(dx/cx)*cx; // ベクトルに周期境界条件適用
			  	dy = dy-rint(dy/cy)*cy;
			  	dz = dz-rint(dz/cz)*cz;
			  
				r = sqrt(dx*dx+dy*dy+dz*dz);    // 原子ijの距離[Å]
				
				// 異種原子間のLJパラメータの計算
				sigIJ = (sig[sp[i]]+sig[sp[j]])*0.5; 
				epsIJ = sqrt(eps[sp[i]]*eps[sp[j]]); 
				s = sigIJ/r;                			
					
				if(r<cutoff){   // r<(カットオフ)の時のみ力を計算
				
					
	             F = 4.0*epsIJ*(1./sigIJ)*   // LJ力の計算 [J/Å]
				           (-12.0*pow(s,13.)+6.0*pow(s,7.));
					                           
				 fx[i] += F*dx/r;     // 原子iについての力の加算
				 fy[i] += F*dy/r;           
				 fz[i] += F*dz/r;           					
					
				 fx[j] -= F*dx/r;     // 原子jについての力の加算
				 fy[j] -= F*dy/r;           
				 fz[j] -= F*dz/r;         					
				}
				
				// ポテンシャルエネルギー[J/mol]の計算
				PE += 4.0*epsIJ*(pow(s,12.)-pow(s,6.));  
				
			}                      // jについてのfor文を閉じる
		}                          // iについてのfor文を閉じる
						
//=======(�F 速度と位置の更新）=======================================

		for(i=0;i<N;i++){           // 原子iについてのfor文
		
         // 速度の更新、単位は[m/s]
			vx[i] += (dt*1.e-15)*(fx[i]*1.e+10)/(w[i]/AVO*0.001); 
			vy[i] += (dt*1.e-15)*(fy[i]*1.e+10)/(w[i]/AVO*0.001); 
			vz[i] += (dt*1.e-15)*(fz[i]*1.e+10)/(w[i]/AVO*0.001); 
	
		 // 座標の更新、単位は[Å]
			rx[i] += vx[i]*(dt*1.e-15)*1.e+10;     
			ry[i] += vy[i]*(dt*1.e-15)*1.e+10; 
			rz[i] += vz[i]*(dt*1.e-15)*1.e+10; 
			
		 // MSD座標の更新、単位は[Å]
			ax[i] += vx[i]*(dt*1.e-15)*1.e+10;     
			ay[i] += vy[i]*(dt*1.e-15)*1.e+10; 
			az[i] += vz[i]*(dt*1.e-15)*1.e+10; 
			
//=======(�G 三次元周期境界条件の適用）==============================

			rx[i] = rx[i] - rint(rx[i]/cx-0.5)*cx;
			ry[i] = ry[i] - rint(ry[i]/cy-0.5)*cy;
			rz[i] = rz[i] - rint(rz[i]/cz-0.5)*cz;

		}  // �Fのiについてのfor文を閉じる


//=======(�H 温度の計算）============================================
								
		KE = 0.0;
		for(i=0;i<N;i++){            
			KE += 
		 (w[i]/AVO*0.001)*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
		}
		
		T = 0.5*KE/(1.5*KB*N); // 温度T [K]
		KE *= 0.5*AVO/N;       // KE:運動エネルギー [J/mol]
		
//=======(�I 温度の調整と重心補正）==================================

		if(k%10==0){             // 温度の調整は10step 毎に行う
			r = sqrt(Tset/T);    // スケーリング係数を計算			
			for(i=0;i<N;i++){
			vx[i] *= r;
			vy[i] *= r;
			vz[i] *= r;
			}
		}
												
		if(k%10==0){                 // 重心補正は10step 毎に行う
			x=y=z=s=0.0;             // 変数の初期化
			for(i=0;i<N;i++){        // 運動量成分(x,y,z)を計算
				x += vx[i]*w[i];  
				y += vy[i]*w[i]; 
				z += vz[i]*w[i];
				s += w[i]; 
			}
			for(i=0;i<N;i++){     //過剰の運動量を各分子から差し引く
				vx[i] = vx[i]-x/(N*s); 
				vy[i] = vy[i]-y/(N*s); 
				vz[i] = vz[i]-z/(N*s); 
			}
		}

//=======(�J 計算経過の出力）========================================
  
		PE = PE*AVO/N;   // 単位変換→[J/mol]
		
		double msdsum[spnum];
					
		if(k==0){            // アニメーション用の初期データ出力
		  fpos <<"1 "<<N<<" "<<(maxstep/10)<<endl;     
		  fpos <<cx<<" "<<cy<<" "<<cz<<endl;         
		  fpos <<"0 "<<dt*10*0.001<<endl;                   
		}

		if(k%10 == 0){         // 結果の出力は10step 毎に行う

		  cout<<k<<"step, KE:"<<KE<<", PE: "<<PE<<", Tem:"<<T
			                   <<" K+P = "<<KE+PE<<endl;
							   
		  fene<<k*dt<<", "<<KE+PE<<", "<<KE<<", "<<PE<<", "<<T<<endl;

		  for(i=0;i<N;i++){       
			fpos << rx[i] <<" "<< ry[i]<<" "<< rz[i]<< endl;
		  }
			
		  fmsd << k*dt<<",  ";
	      for(j=0;j<spnum;j++){  // 各成分ごとにMSDを計算する
			msdsum[j] = 0.0;
			m         =   0;
			for(i=0;i<N;i++){
			 if(sp[i]==j){
			   msdsum[j] += ax[i]*ax[i]+ay[i]*ay[i]+az[i]*az[i];
				m++;
			 }
			}
			fmsd << msdsum[j]/m <<", ";
		  }
		  fmsd << endl;
		}
		
//=======(�K 反復回数の終了判定）====================================

	}  // kについてのfor文を閉じる
	
//=======(�L 計算終了）==============================================
			
	fpos.close();      // ファイルを閉じる
	fene.close();
	fmsd.close();	
	
	cout << "MD計算が全て終了しました。"<<endl;
			
}                      //main文を閉じる

//===================================================================