这个程序运行停止 高手来看看
// EMD.cpp : Defines the entry point for the console application.//功能：采用平衡态分子动力学算法计算Ar的热传导性能
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <iostream>
//#include <mpi.h>
#define Si 1
#define other 2
//定义全局变量
//MD常数定义
//SW势函数中参数
//硅与硅之间的参数
//数组常数
extern "C" const int nmax=20000, listmax=2000000;
double eps2_De,epsavg_De,sig2_Re,sigavg_Re;
//计算用变量
int rank,p;
int n,ncx,ncy,ncz,nSi;
double tinit; //初始温度
double Mass_tot;
//粒子的位置、速度、受力
double rx0[nmax],ry0[nmax],rz0[nmax];
int typeatom[nmax];
//记录每个粒子的势能
double poten[nmax];
//计算两体有效热流用
double afx[listmax],afy[listmax],afz[listmax];
double arx[listmax],ary[listmax],arz[listmax];
//计算三体有效热流用
int f3count;
int atomi[listmax*3],atomj[listmax*3],atomk[listmax*3];
double afx3[listmax*3],afy3[listmax*3],afz3[listmax*3];
//构件邻居列表用
int pointers[listmax],neighbor[listmax];
//模拟步数等
int timestep,currentstep,totsteps,startmeas,freq;
//模拟区域大小、及其容积
double boxx,boxy,boxz,volume;
//温度、有效热流的均值
double tavg, jxavg,jyavg,jzavg;
//标示相变的参数，接近1时为固相
double topr;
//总能量、动能、势能、温度、动量
double tote,ke,pe,T,px,py,pz;
//有效热流及其分量
double jx,jy,jz,jimx,jimy,jimz,jsdx ,jsdy,jsdz;
//文件指针
FILE *totflux,*temperature,*imsdflx,*cnstot,*potential,*energy,*datafile;
FILE *avgs,*position,*initpara,*velocity,*force;
//=='y' or 'Y'，继续上次计算，否则重新计算
void setfcc()
{
double xbase[]={0.0, 0.5, 0.5, 0.0, 0.25, 0.25, 0.75, 0.75};
double ybase[]={0.0, 0.5, 0.0, 0.5, 0.25, 0.75, 0.25, 0.75};
double zbase[]={0.0, 0.0, 0.5, 0.5, 0.25, 0.75, 0.75, 0.25};
int counter, counter1;
counter=0, counter1=0;
for(int ix=0; ix<ncx; ix++)
{
for(int iy=0; iy<ncy; iy++)
{
for(int iz=0; iz<ncz; iz++)
{
for(int i=0; i<8; i++)
{
counter++;
typeatom[counter]=Si;
rx[counter]=xbase[i]*a+ix*a;
ry[counter]=ybase[i]*a+iy*a;
rz[counter]=zbase[i]*a+iz*a;
}
}
}
}
typeatom[336]=other;
double RRx,RRy,RRz,k;
int ii=0;
for(int i=0;i<n;i++)
{
RRx=rx[i]-rx[336];
RRy=ry[i]-ry[336];
RRz=rz[i]-rz[336];
// if(RRx<0 || RRy<0 || RRz<0) continue;
// if(RRx==0 && RRy==0 && RRz==0) continue;
if(fabs(RRx)<0.5&&fabs(RRy)<0.5&&fabs(RRz)<0.5)
{
typeatom[i]=other;
counter1++;
}
if(counter1>=19) break;
}
for(int jj=0;jj<n;jj++)
{
if(typeatom[jj]==other) ii++;
}
nSi=n-ii;
printf("n=%d\n",n); //粒子总数
printf("ii=%d\n",ii); //杂质数量
printf("counter=%d\n",counter);
printf("counter1=%d\n",counter1);
k=ii/512.0; //杂质浓度
printf("k=%lf\n",k);
if(counter!=(n))
error("setfcc: counter/=n",0);
}
double uniform()
{
int L, C, M;
static int SEED=0;
L = 1029;
C = 221591;
M = 1048576;
SEED = (SEED * L+C)%M;
return double(SEED)/M;
}
double gauss()
{
double a1,a3,a5,a7,a9;
a1=3.949846138;
a3 = 0.252408784;
a5 = 0.076542912;
a7 = 0.008355968;
a9 = 0.029899776;
double sum=0.0;
for(int i=0;i<12; i++)
sum=sum+uniform();
double r,r2,res;
r=(sum-6.0)/4.0;
r2=r*r;
res = ((((a9*r2+a7)*r2+a5)*r2+a3)*r2+a1)*r;
return res;
}
void initvel()
{
double pxt,pyt,pzt,factor;
pxt=0.0;
pyt=0.0;
pzt=0.0;
int i;
factor=sqrt(tinit);
for( i=0; i<n; i++)
{
if(typeatom[i]==Si)
{
vx[i]=factor*gauss();
vy[i]=factor*gauss();
vz[i]=factor*gauss();
pxt=pxt+vx[i];
pyt=pyt+vy[i];
pzt=pzt+vz[i];
}
else
{
vx[i]=factor*gauss()/sqrt(mass21);
vy[i]=factor*gauss()/sqrt(mass21);
vz[i]=factor*gauss()/sqrt(mass21);
//净动量
pxt=pxt+vx[i]*mass21;
pyt=pyt+vy[i]*mass21;
pzt=pzt+vz[i]*mass21;
}
}
Mass_tot=(n-nSi)*mass21+nSi;
double VsubX=pxt/Mass_tot;
double VsubY=pyt/Mass_tot;
double VsubZ=pzt/Mass_tot;
//Velsq=∑(m*v*v)，即2倍的动能
double VelsqX=0.0,VelsqY=0.0,VelsqZ=0.0;
for(i=0; i<n; i++)
{
//将系统的总动量置零
vx[i]=vx[i]-VsubX;
vy[i]=vy[i]-VsubY;
vz[i]=vz[i]-VsubZ;
if(typeatom[i]==Si)
{
VelsqX=VelsqX+vx[i]*vx[i];
VelsqY=VelsqY+vy[i]*vy[i];
VelsqZ=VelsqZ+vz[i]*vz[i];
}
else
{
VelsqX=VelsqX+vx[i]*vx[i]*mass21;
VelsqY=VelsqY+vy[i]*vy[i]*mass21;
VelsqZ=VelsqZ+vz[i]*vz[i]*mass21;
}
}
//保证X,Y,Z三个方向的动能相同
//由于动量是自动守恒的，故这里的自由度为n-1.
double factorx,factory,factorz;
factorx=sqrt((n-1)*tinit/VelsqX);
factory=sqrt((n-1)*tinit/VelsqY);
factorz=sqrt((n-1)*tinit/VelsqZ);
for(i=0; i<n; i++)
{
vx[i]=vx[i]*factorx;
vy[i]=vy[i]*factory;
vz[i]=vz[i]*factorz;
}
}
void applypbc(double &rxij, double &ryij, double &rzij)
{
if(rxij>boxx*0.5)
rxij=rxij-boxx;
if(rxij<-boxx*0.5)
rxij=rxij+boxx;
if(ryij>boxy*0.5)
ryij=ryij-boxy;
if(ryij<-boxy*0.5)
ryij=ryij+boxy;
if(rzij>boxz*0.5)
rzij=rzij-boxz;
if(rzij<-boxz*0.5)
rzij=rzij+boxz;
}
void save()
{
for(int i=0; i<n; i++)
{
rx0[i]=rx[i];
ry0[i]=ry[i];
rz0[i]=rz[i];
}
}
void update()
{
double rxi,ryi,rzi,rxij,ryij,rzij,rijsq;
double rlistsq;
rlistsq=rlist*rlist;
int pair=0;
if(pair>listmax)
error("update: listmax is too small\n",0);
for(int i=0; i<n; i++)
{
// rx0[i]=rx[i];
// ry0[i]=ry[i];
// rz0[i]=rz[i];
pointers[i]=pair;
rxi=rx[i];
ryi=ry[i];
rzi=rz[i];
for(int j=0; j<n; j++)
{
rxij=rxi-rx[j];
ryij=ryi-ry[j];
rzij=rzi-rz[j];
//pbc
if(rxij>boxx*0.5)
rxij=rxij-boxx;
if(rxij<-boxx*0.5)
rxij=rxij+boxx;
if(ryij>boxy*0.5)
ryij=ryij-boxy;
if(ryij<-boxy*0.5)
ryij=ryij+boxy;
if(rzij>boxz*0.5)
rzij=rzij-boxz;
if(rzij<-boxz*0.5)
rzij=rzij+boxz;
rijsq=rxij*rxij+ryij*ryij+rzij*rzij;
if(rijsq<rlistsq)
{
neighbor[pair]=j;
pair++;
if(pair>listmax)
error("update: listmax is too small\n",0);
}
}
}
pointers[n]=pair;
}
bool checkcomp()
{
double dispmax;
dispmax=0.0;
for(int i=0; i<n; i++)
{
if(fabs(rx[i]-rx0[i])>dispmax) dispmax=fabs(rx[i]-rx0[i]);
if(fabs(ry[i]-ry0[i])>dispmax) dispmax=fabs(ry[i]-ry0[i]);
if(fabs(rz[i]-rz0[i])>dispmax) dispmax=fabs(rz[i]-rz0[i]);
}
dispmax=2.0*sqrt(dispmax*dispmax*3.0);
if(dispmax>(rlist-rcut))
return true;
else
return false;
}
void swforce()
{
//内存清零
memset(fx,0,sizeof(double)*nmax);
memset(fy,0,sizeof(double)*nmax);
memset(fz,0,sizeof(double)*nmax);
memset(poten,0,sizeof(double)*nmax);
pe=0.0;
double rxij,ryij,rzij,rijsq,rij;
double fxij,fyij,fzij;
double fij,potenij;
double rcutsq=rcut*rcut-0.001;
int i,j,k;
//计算两体作用力
for(i=0; i<n-1; i++)
{
for(k=pointers[i]; k<pointers[i+1]; k++)
{
j=neighbor[k];
afx[k]=0.0;
afy[k]=0.0;
afz[k]=0.0;
//只计算与编号大于自身的邻居之间的作用力
if(i>j) continue;
rxij=rx[i]-rx[j];
ryij=ry[i]-ry[j];
rzij=rz[i]-rz[j];
applypbc(rxij,ryij,rzij);
rijsq=rxij*rxij+ryij*ryij+rzij*rzij;
if(rijsq>rcutsq) continue;
rij=sqrt(rijsq);
if((typeatom[i]==1)&&( typeatom [j]==1)) //两粒子同为Si原子
{
fij=par_A_si/rij*(4.0*par_B_si/pow(rij,5)+(par_B_si/pow(rij,4)-1.0)/pow(rij-rcut,2))*exp(1.0/(rij-rcut));
//计算力在各方向上的分量
fxij=fij*rxij;
fyij=fij*ryij;
fzij=fij*rzij;
fx[i]=fx[i]+fxij;
fy[i]=fy[i]+fyij;
fz[i]=fz[i]+fzij;
potenij=par_A_si*(par_B_si/pow(rij,4)-1.0)*exp(1.0/(rij-rcut));
poten[i]=poten[i]+potenij*0.5;
pe=pe+potenij;
fx[j]=fx[j]-fxij;
fy[j]=fy[j]-fyij;
fz[j]=fz[j]-fzij;
poten[j]=poten[j]+potenij*0.5;
afx[k]=fxij;
afy[k]=fyij;
afz[k]=fzij;
}
else if((typeatom[i]==2)&&( typeatom [j]==2)) //两粒子同为Ge原子 1.04=sig_ge/sig_si
{
fij=par_A_ge/rij*(4.0*par_B_ge*1.04/pow(rij,5)+(par_B_ge/pow(rij,4)-1.0)*pow(1.04,2)/pow(rij-rcut,2))*exp(1.04/(rij-rcut));
//计算力在各方向上的分量
fxij=fij*rxij;
fyij=fij*ryij;
fzij=fij*rzij;
fx[i]=fx[i]+fxij;
fy[i]=fy[i]+fyij;
fz[i]=fz[i]+fzij;
potenij=par_A_ge*(par_B_ge/pow(rij,4)-1.0)*exp(1.04/(rij-rcut));
poten[i]=poten[i]+potenij*0.5;
pe=pe+potenij;
fx[j]=fx[j]-fxij;
fy[j]=fy[j]-fyij;
fz[j]=fz[j]-fzij;
poten[j]=poten[j]+potenij*0.5;
afx[k]=fxij;
afy[k]=fyij;
afz[k]=fzij;
}
else if(((typeatom[i]==1)&&( typeatom [j]==2))||((typeatom[i]==2)&&( typeatom [j]==1))) //Si Ge原子间相互作用!! 1.04=sig_avg/sig_si
{
fij=par_A_avg/rij*(4.0*par_B_avg*1.02/pow(rij,5)+(par_B_avg/pow(rij,4)-1.0)*pow(1.02,2)/pow(rij-rcut,2))*exp(1.02/(rij-rcut));
//计算力在各方向上的分量
fxij=fij*rxij;
fyij=fij*ryij;
fzij=fij*rzij;
fx[i]=fx[i]+fxij;
fy[i]=fy[i]+fyij;
fz[i]=fz[i]+fzij;
potenij=par_A_avg*(par_B_avg/pow(rij,4)-1.0)*exp(1.02/(rij-rcut));
poten[i]=poten[i]+potenij*0.5;
pe=pe+potenij;
fx[j]=fx[j]-fxij;
fy[j]=fy[j]-fyij;
fz[j]=fz[j]-fzij;
poten[j]=poten[j]+potenij*0.5;
afx[k]=fxij;
afy[k]=fyij;
afz[k]=fzij;
}
}
//计算三体作用力
double dxij,dyij,dzij,dxik,dyik,dzik;
double rxik,ryik,rzik,riksq,rik;
double fcvaljx,fcvalkx,fcvaljy,fcvalky,fcvaljz,fcvalkz;
double cr,er, pij,pik;
double potenijk;
f3count=-1;
for(i=0; i<n; i++)
{
//只计算H(i,j,k),对于U3=H(ijk)+H(jki)+H(kij)分别在三次循环中计算
for(int second=pointers[i]; second<pointers[i+1]-1; second++)
{
j=neighbor[second];
rxij=rx[i]-rx[j];
ryij=ry[i]-ry[j];
rzij=rz[i]-rz[j];
applypbc(rxij,ryij,rzij);
rijsq=rxij*rxij+ryij*ryij+rzij*rzij;
if(rijsq>rcutsq) continue;
rij=sqrt(rijsq);
dxij=rxij/rij;
dyij=ryij/rij;
dzij=rzij/rij;
for(int third=second+1; third<pointers[i+1]; third++)
{
k=neighbor[third];
rxik=rx[i]-rx[k];
ryik=ry[i]-ry[k];
rzik=rz[i]-rz[k];
applypbc(rxik,ryik,rzik);
riksq=rxik*rxik+ryik*ryik+rzik*rzik;
if(riksq>rcutsq) continue;
rik=sqrt(riksq);
dxik=rxik/rik;
dyik=ryik/rik;
dzik=rzik/rik;
cr=dxij*dxik+dyij*dyik+dzij*dzik; //cosθ213
if((typeatom[i]==1)&&(typeatom[j]==1)&&(typeatom[k]==1))//三种粒子同是Si
{
er=par_lamda_si*(cr+1.0/3.0)*exp(par_gama_si/(rij-rcut)+par_gama_si/(rik-rcut));
pij=par_gama_si*(cr+1.0/3.0)/pow(rij-rcut,2);
pik=par_gama_si*(cr+1.0/3.0)/pow(rik-rcut,2);
fcvaljx=-er*(pij*dxij+2.0*(cr*dxij-dxik)/rij);
fcvalkx=-er*(pik*dxik+2.0*(cr*dxik-dxij)/rik);
fx[i]= fx[i]-fcvaljx-fcvalkx;
fx[j]= fx[j]+fcvaljx;
fx[k]=fx[k]+fcvalkx;
fcvaljy=-er*(pij*dyij+2.0*(cr*dyij-dyik)/rij);
fcvalky=-er*(pik*dyik+2.0*(cr*dyik-dyij)/rik);
fy[i]= fy[i]-fcvaljy-fcvalky;
fy[j]= fy[j]+fcvaljy;
fy[k]= fy[k]+fcvalky;
fcvaljz=-er*(pij*dzij+2.0*(cr*dzij-dzik)/rij);
fcvalkz=-er*(pik*dzik+2.0*(cr*dzik-dzij)/rik);
fz[i]= fz[i]-fcvaljz-fcvalkz;
fz[j]= fz[j]+fcvaljz;
fz[k]= fz[k]+fcvalkz;
potenijk=er*(cr+1.0/3.0);
//根据能量均分理论
poten[i]= poten[i]+potenijk/3.0;
poten[j]= poten[j]+potenijk/3.0;
poten[k]= poten[k]+potenijk/3.0;
//记录三体作用力对有效热流的贡献
f3count++;
atomi[f3count]=i;//i,j,k均为全局编号
atomj[f3count]=j;
atomk[f3count]=k;
afx3[f3count]=-fcvaljx-fcvalkx;
afy3[f3count]=-fcvaljy-fcvalky;
afz3[f3count]=-fcvaljz-fcvalkz;
f3count++;
atomi[f3count]=j;
atomj[f3count]=i;
atomk[f3count]=k;
afx3[f3count]=fcvaljx;
afy3[f3count]=fcvaljy;
afz3[f3count]=fcvaljz;
f3count++;
atomi[f3count]=k;
atomj[f3count]=i;
atomk[f3count]=j;
afx3[f3count]=fcvalkx;
afy3[f3count]=fcvalky;
afz3[f3count]=fcvalkz;
}
else if((typeatom[i]==2)&&(typeatom[j]==2)&&(typeatom[k]==2))//三个粒子同为Ge
{
er=par_lamda_ge*(cr+1.0/3.0)*exp(par_gama_ge/(rij-rcut)+par_gama_ge/(rik-rcut));
pij=par_gama_ge*1.04*(cr+1.0/3.0)/pow(rij-rcut,2);
pik=par_gama_ge*1.04*(cr+1.0/3.0)/pow(rik-rcut,2);
fcvaljx=-er*(pij*dxij+2.0*1.04*(cr*dxij-dxik)/rij);
fcvalkx=-er*(pik*dxik+2.0*1.04*(cr*dxik-dxij)/rik);
fx[i]= fx[i]-fcvaljx-fcvalkx;
fx[j]= fx[j]+fcvaljx;
fx[k]= fx[k]+fcvalkx;
fcvaljy=-er*(pij*dyij+2.0*1.04*(cr*dyij-dyik)/rij);
fcvalky=-er*(pik*dyik+2.0*1.04*(cr*dyik-dyij)/rik);
fy[i]= fy[i]-fcvaljy-fcvalky;
fy[j]= fy[j]+fcvaljy;
fy[k]= fy[k]+fcvalky;
fcvaljz=-er*(pij*dzij+2.0*1.04*(cr*dzij-dzik)/rij);
fcvalkz=-er*(pik*dzik+2.0*1.04*(cr*dzik-dzij)/rik);
fz[i]= fz[i]-fcvaljz-fcvalkz;
fz[j]= fz[j]+fcvaljz;
fz[k]= fz[k]+fcvalkz;
potenijk=er*(cr+1.0/3.0);
poten[i]= poten[i]+potenijk/3.0;
poten[j]= poten[j]+potenijk/3.0;
poten[k]= poten[k]+potenijk/3.0;
//记录三体作用力对有效热流的贡献
f3count++;
atomi[f3count]=i;//i,j,k均为全局编号
atomj[f3count]=j;
atomk[f3count]=k;
afx3[f3count]=-fcvaljx-fcvalkx;
afy3[f3count]=-fcvaljy-fcvalky;
afz3[f3count]=-fcvaljz-fcvalkz;
f3count++;
atomi[f3count]=j;
atomj[f3count]=i;
atomk[f3count]=k;
afx3[f3count]=fcvaljx;
afy3[f3count]=fcvaljy;
afz3[f3count]=fcvaljz;
f3count++;
atomi[f3count]=k;
atomj[f3count]=i;
atomk[f3count]=j;
afx3[f3count]=fcvalkx;
afy3[f3count]=fcvalky;
afz3[f3count]=fcvalkz;
}
else
{
er=par_lamda_avg*(cr+1.0/3.0)*exp(par_gama_avg/(rij-rcut)+par_gama_avg/(rik-rcut));
pij=par_gama_avg*1.02*(cr+1.0/3.0)/pow(rij-rcut,2);
pik=par_gama_avg*1.02*(cr+1.0/3.0)/pow(rik-rcut,2);
fcvaljx=-er*(pij*dxij+2.0*1.02*(cr*dxij-dxik)/rij);
fcvalkx=-er*(pik*dxik+2.0*1.02*(cr*dxik-dxij)/rik);
fx[i]= fx[i]-fcvaljx-fcvalkx;
fx[j]= fx[j]+fcvaljx;
fx[k]= fx[k]+fcvalkx;
fcvaljy=-er*(pij*dyij+2.0*1.02*(cr*dyij-dyik)/rij);
fcvalky=-er*(pik*dyik+2.0*1.02*(cr*dyik-dyij)/rik);
fy[i]= fy[i]-fcvaljy-fcvalky;
fy[j]= fy[j]+fcvaljy;
fy[k]= fy[k]+fcvalky;
fcvaljz=-er*(pij*dzij+2.0*1.02*(cr*dzij-dzik)/rij);
fcvalkz=-er*(pik*dzik+2.0*1.02*(cr*dzik-dzij)/rik);
fz[i]= fz[i]-fcvaljz-fcvalkz;
fz[j]= fz[j]+fcvaljz;
fz[k]= fz[k]+fcvalkz;
potenijk=er*(cr+1.0/3.0);
poten[i]= poten[i]+potenijk/3.0;
poten[j]= poten[j]+potenijk/3.0;
poten[k]= poten[k]+potenijk/3.0;
//记录三体作用力对有效热流的贡献
f3count++;
atomi[f3count]=i;//i,j,k均为全局编号
atomj[f3count]=j;
atomk[f3count]=k;
afx3[f3count]=-fcvaljx-fcvalkx;
afy3[f3count]=-fcvaljy-fcvalky;
afz3[f3count]=-fcvaljz-fcvalkz;
f3count++;
atomi[f3count]=j;
atomj[f3count]=i;
atomk[f3count]=k;
afx3[f3count]=fcvaljx;
afy3[f3count]=fcvaljy;
afz3[f3count]=fcvaljz;
f3count++;
atomi[f3count]=k;
atomj[f3count]=i;
atomk[f3count]=j;
afx3[f3count]=fcvalkx;
afy3[f3count]=fcvalky;
afz3[f3count]=fcvalkz;
}
}
}
}
}
}
//采用velocity verlet积分算法，对位移、速度进行积分。
//位移前进一步，速度前进半步
void move1()
{
double dtsq;
dtsq=dt*dt;
for(int i=0; i<n; i++)
{
if(typeatom[i]==Si)
{
rx[i]=rx[i]+dt*vx[i]+dtsq*fx[i]*0.5;
ry[i]=ry[i]+dt*vy[i]+dtsq*fy[i]*0.5;
rz[i]=rz[i]+dt*vz[i]+dtsq*fz[i]*0.5;
vx[i]=vx[i]+fx[i]*dt*0.5;
vy[i]=vy[i]+fy[i]*dt*0.5;
vz[i]=vz[i]+fz[i]*dt*0.5;
}
else if(typeatom[i]==other)
{
rx[i]=rx[i]+dt*vx[i]+dtsq*fx[i]*0.5/mass21;
ry[i]=ry[i]+dt*vy[i]+dtsq*fy[i]*0.5/mass21;
rz[i]=rz[i]+dt*vz[i]+dtsq*fz[i]*0.5/mass21;
vx[i]=vx[i]+fx[i]*dt*0.5/mass21;
vy[i]=vy[i]+fy[i]*dt*0.5/mass21;
vz[i]=vz[i]+fz[i]*dt*0.5/mass21;
}
}
}
//速度前进半步
void move2()
{
for(int i=0; i<n; i++)
{
if(typeatom[i]==Si)
{
vx[i]=vx[i]+fx[i]*dt*0.5;
vy[i]=vy[i]+fy[i]*dt*0.5;
vz[i]=vz[i]+fz[i]*dt*0.5;
}
else if(typeatom[i]==other)
{
vx[i]=vx[i]+fx[i]*dt*0.5/mass21;
vy[i]=vy[i]+fy[i]*dt*0.5/mass21;
vz[i]=vz[i]+fz[i]*dt*0.5/mass21;
}
}
}
void sumkemom()
{
px=0.0;
py=0.0;
pz=0.0;
ke=0.0;
T=0.0;
for(int i=0; i<n; i++)
{
if(typeatom[i]==Si)
{
px=px+vx[i];
py=py+vy[i];
pz=pz+vz[i];
ke=ke+vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
}
else if(typeatom[i]==other)
{
px=px+mass21*vx[i];
py=py+mass21*vy[i];
pz=pz+mass21*vz[i];
ke=ke+mass21*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
}
}
T=ke/(3.0*n);
ke=ke*0.5;
if(timestep>startmeas)
tavg=tavg+T;
}
//注意这里使用到系统的动能、动量。因此，系统的动能、动量必须在此之前求出。
void adjusttemperature()
{
double Psq=px*px+py*py+pz*pz;
//临时变量
double TempA=2.0*ke-Psq/Mass_tot;
double TempB=3.0*n*tinit-Psq/Mass_tot;
double factor=sqrt(TempB/TempA);
double VsubX=(factor-1.0)*px/Mass_tot;
double VsubY=(factor-1.0)*py/Mass_tot;
double VsubZ=(factor-1.0)*pz/Mass_tot;
//对于所有原子进行速度的调整
//？？？但是对于不同质量的原子进行同比例的速度的缩放是否合适？？？
for(int i=0; i<n; i++)
{
vx[i]=vx[i]*factor-VsubX;
vy[i]=vy[i]*factor-VsubY;
vz[i]=vz[i]*factor-VsubZ;
}
}
void getflux()
{
jimx=0.0;
jimy=0.0;
jimz=0.0;
jsdx=0.0;
jsdy=0.0;
jsdz=0.0;
int pairstart,pairstop,j;
double fdotv,velsqi;
for(int i=0; i<n; i++)
{
pairstart=pointers[i];
pairstop=pointers[i+1]-1;
if(pairstart<=pairstop)
{
for(int k=pairstart; k<=pairstop;k++)
{
j=neighbor[k];
fdotv=afx[k]*vx[i]+afy[k]*vy[i]+afz[k]*vz[i];
jimx=jimx+arx[k]*fdotv;
jimy=jimy+ary[k]*fdotv;
jimz=jimz+arz[k]*fdotv;
fdotv=afx[k]*vx[j]+afy[k]*vy[j]+afz[k]*vz[j];
jimx=jimx+arx[k]*fdotv;
jimy=jimy+ary[k]*fdotv;
jimz=jimz+arz[k]*fdotv;
}
}
if(typeatom[i]==Si)
velsqi=(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i])*0.5;
else if(typeatom[i]==other)
velsqi=mass21*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i])*0.5;
jsdx=jsdx+(velsqi+poten[i])*vx[i];
jsdy=jsdy+(velsqi+poten[i])*vy[i];
jsdz=jsdz+(velsqi+poten[i])*vz[i];
}
jimx=jimx*0.5/volume;
jimy=jimy*0.5/volume;
jimz=jimz*0.5/volume;
jsdx=jsdx/volume;
jsdy=jsdy/volume;
jsdz=jsdz/volume;
jx=jimx+jsdx;
jy=jimy+jsdy;
jz=jimz+jsdz;
jxavg=jxavg+jx;
jyavg=jyavg+jy;
jzavg=jzavg+jz;
}
int main(int argc, char* argv[])
{
// MPI_Init(&argc,&argv); //MPI程序初始化
// MPI_Comm_rank(MPI_COMM_WORLD,&rank);
// MPI_Comm_size(MPI_COMM_WORLD,&p);
openfiles();
readinitpara();
if(totsteps<=startmeas)
error("totsteps<=startmeas", 0);
//温度无量纲化
tinit=tinit*kb/eps;
eps2_De=eps2/eps;
epsavg_De=epsavg/eps;
sig2_Re=sig2/sig;
sigavg_Re=sigavg/sig;
n=8*ncx*ncy*ncz; //总粒子数
// cout<<"粒子总数为:"<<n<<endl;
if(nmax<n)
error("nmax<n", 0);
//模拟区域的大小
boxx=a*ncx;
boxy=a*ncy;
boxz=a*ncz;
volume=boxx*boxy*boxz;
writeinitpara();
//初始化参数
tavg=0.0;
jxavg=0.0;
jyavg=0.0;
jzavg=0.0;
timestep=1;
currentstep=1;
setfcc();
writepos();
initvel();
writevel();
save();
update();
swforce();
writeforce();
writepotential();
sumkemom();
tote=ke+pe;
writetemperature();
writeenergy();
//时间演进
bool res;
for(timestep=currentstep+1; timestep<=totsteps;timestep++)
{
if(timestep%1000==0) printf("%7d\n",timestep);
move1();
res=checkcomp();
if(res)
{
save();
update();
}
swforce();
move2();
sumkemom();
tote=ke+pe;
if(timestep%freq==0)
{
writeenergy();
writetemperature();
writepotential();
}
//注意：sumkemom必须在adjusttemperature之前执行
if((timestep<startmeas/2)&&(timestep%5==0))
adjusttemperature();
if(timestep%50000==0)
{
writepos();
writevel();
writeforce();
}
//从startmeas开始，每步计算一次有效热流
if(timestep>=startmeas)
{
getflux();
writeflux();
}
/*
if(timestep%100==0)
{
top();
writetop();
}
*/
}
if(currentstep>startmeas)
startmeas=currentstep;
tavg=tavg/(totsteps-startmeas);
jxavg=jxavg/(totsteps-startmeas);
jyavg=jyavg/(totsteps-startmeas);
jzavg=jzavg/(totsteps-startmeas);
writeavgs();
return 0;
// MPI_Finalize();
}
