汉化的MD小程序

/******************************************************************************* 兰纳-琼斯势： 又称L-J势能函数或6-12势能函数，用来描述惰性气体分子间相互作用， 近距离排斥，远距离吸引，本身是个近似公式。 V（r）=e*((r_min/r)^12-(r_min/r)^6)),其中e是势能阱的深度，r_min是相互作用势能 正好为0时的距离。 用法 %md < md.in (输入文件md.in的格式见md.h) *******************************************************************************/ #include <stdio.h> #include <math.h> #include "md.h" #include <time.h> #include <iostream> #include<fstream> using namespace std; main(int argc, char **argv) { double cpu, cpu1, cpu2; fstream fp,outfile; char filename[50]; double vSum,vv; int n,k; if (argc<2) { cout<<"请输入文件名："<<endl; cin>>filename;} cout<<"argv="<<filename<<endl; fp.open(filename, ios::in); if(!fp) cout<<"无法打开文件"<<filename<<endl; /* 读入文件in的8个数据 */ fp>> InitUcell[0]>>InitUcell[1]>>InitUcell[2]; fp>> Density >>InitTemp>>DeltaT>>StepLimit>>StepAvg; fp.close(); /*--------以上读入文件in的数据-------------------------------------------*/ outfile.open("md.out",ios::out); InitParams(); InitConf(); ComputeAccel(); /* 计算初始加速度 */ cpu1 = ((double) clock())/CLOCKS_PER_SEC; /* 计算时间单位秒 */ /* clock函数指程序启动至调用本函数时的计时单元数，为long 型， 而一秒内含有多少计时单元用常数 CLOCKS_PER_SEC */ for (stepCount=1; stepCount<=StepLimit; stepCount++) { SingleStep(); ApplyBoundaryCond(); if (stepCount%StepAvg == 0) { /*------------------------------------------------------------------------------ 估算物理参数: 动能、势能及总能. ------------------------------------------------------------------------------*/ vSum = kinEnergy = 0.0; for (n=0; n<nAtom; n++) { vv = 0.0; for (k=0; k<3; k++) { vSum = vSum + rv[n][k]; /* 三维速度和 */ vv = vv + rv[n][k]*rv[n][k]; /* 三维速度平方和 */ } kinEnergy = kinEnergy + vv; /* 所有原子三维速度平方和累加 */ } kinEnergy *= (0.5/nAtom); /* 除以2，除以原子数 */ potEnergy /= nAtom; /* 平均势能 */ totEnergy = kinEnergy + potEnergy; /* 总能=动能+势能 */ temperature = kinEnergy*2.0/3.0; /* 温度= 动能/3*2 */ /* 输出上述参数 */ outfile<<stepCount*DeltaT<<" "<<temperature<<" "<<potEnergy<<" "<<totEnergy<<endl; } } cpu2 = ((double) clock())/CLOCKS_PER_SEC; cpu = cpu2-cpu1; outfile<<"cpu1="<<cpu1<<" cpu2="<<cpu2<<" cpu="<< cpu<<endl; outfile.close(); return 0; } /*----------------------------------------------------------------------------*/ void InitParams() { /*------------------------------------------------------------------------------ 初始化参数 ------------------------------------------------------------------------------*/ int k; double rr,ri2,ri6,r1; /* 计算基础数据，类型均在H文件中声明过 */ DeltaTH = 0.5*DeltaT; /* DeltaTH是半步长 */ for (k=0; k<3; k++) { Region[k] = InitUcell[k]/pow(Density/4.0,1.0/3.0); /* 盒子三维长度计算。 实在看不懂！ pow是指数函数，10^4为pow(10,4); */ RegionH[k] = 0.5*Region[k]; /* 半模型长度 */ } /* Constants for potential truncation */ rr = RCUT*RCUT; ri2 = 1.0/rr; ri6 = ri2*ri2*ri2; r1=sqrt(rr); Uc = 4.0*ri6*(ri6 - 1.0); Duc = -48.0*ri6*(ri6 - 0.5)/r1; /* rr是r的平方，ri2是平方倒数，ri6是倒数6次方，r1就是RCUT */ } /*----------------------------------------------------------------------------*/ void InitConf() { /*------------------------------------------------------------------------------ r 初始化为面心立方位置. rv 初始化成与温度相关的随机速度. FCC，即面心立方晶格，是晶体结构的一种。面心立方晶格的晶胞是一个立方体， 立方体的八个顶角和六个面的中心各有一个原子。 ------------------------------------------------------------------------------*/ double c[3],gap[3],e[3],vSum[3],vMag; int j,n,k,nX,nY,nZ; double seed; /* 设置面心立方的原子位置模板 */ double origAtom[4][3] = {{0.0, 0.0, 0.0}, {0.0, 0.5, 0.5}, {0.5, 0.0, 0.5}, {0.5, 0.5, 0.0}}; /* 设置FCC晶格 */ for (k=0; k<3; k++) gap[k] = Region[k]/InitUcell[k]; nAtom = 0; for (nZ=0; nZ<InitUcell[2]; nZ++) { c[2] = nZ*gap[2]; for (nY=0; nY<InitUcell[1]; nY++) { c[1] = nY*gap[1]; for (nX=0; nX<InitUcell[0]; nX++) { c[0] = nX*gap[0]; /* 设置FCC晶格 */ for (j=0; j<4; j++) { for (k=0; k<3; k++) r[nAtom][k] = c[k] + gap[k]*origAtom[j][k]; ++nAtom; /* 此步中c[k]是间隔，gap[k]是步长，origAtom[][]决定方向。 为所有原子设定 晶格位置。 */ } } } } /* 生成随机的速度 */ seed = 13597.0; vMag = sqrt(3*InitTemp); for(k=0; k<3; k++) vSum[k] = 0.0; for(n=0; n<nAtom; n++) { RandVec3(e,&seed); /* 各原子速度随便定 */ for (k=0; k<3; k++) { rv[n][k] = vMag*e[k]; /* 速度跟温度扯蛋？ */ vSum[k] = vSum[k] + rv[n][k]; /* 计算总动量，为归零准备 */ } } /* 保持总动量为0 */ for (k=0; k<3; k++) vSum[k] = vSum[k]/nAtom; for (n=0; n<nAtom; n++) for(k=0; k<3; k++) rv[n][k] = rv[n][k] - vSum[k]; /* 各原子在各个方向上的速度均减去总平均速度 */ } /*----------------------------------------------------------------------------*/ void ComputeAccel() { /*------------------------------------------------------------------------------ ra,加速度，根据原子坐标r,使用Lennard-Jones势能公式计算. ------------------------------------------------------------------------------*/ double dr[3],f,fcVal,rrCut,rr,ri2,ri6,r1; int j1,j2,n,k; /* 首先初始为零 */ rrCut = RCUT*RCUT; for (n=0; n<nAtom; n++) for (k=0; k<3; k++) ra[n][k] = 0.0; potEnergy = 0.0; /* 原子对双嵌套循环 */ for (j1=0; j1<nAtom-1; j1++) { for (j2=j1+1; j2<nAtom; j2++) { /* 计算原子三维上原子距离平方和 */ for (rr=0.0, k=0; k<3; k++) { dr[k] = r[j1][k] - r[j2][k]; /*计算各原子间的距离 */ /* 判断距离与半晶格的差，设成最近的原子 */ dr[k] = dr[k] - SignR(RegionH[k],dr[k]-RegionH[k]) - SignR(RegionH[k],dr[k]+RegionH[k]); rr = rr + dr[k]*dr[k]; } /* 计算极限距离内的加速度和势能 */ if (rr < rrCut) { ri2 = 1.0/rr; ri6 = ri2*ri2*ri2; r1 = sqrt(rr); fcVal = 48.0*ri2*ri6*(ri6-0.5) + Duc/r1; for (k=0; k<3; k++) { /* 计算三向上的加速度 */ f = fcVal*dr[k]; /* 求力 */ ra[j1][k] = ra[j1][k] + f; ra[j2][k] = ra[j2][k] - f; } /* 计 算 动 能 */ potEnergy = potEnergy + 4.0*ri6*(ri6-1.0) - Uc - Duc*(r1-RCUT); } } } } /*----------------------------------------------------------------------------*/ void SingleStep() { /*------------------------------------------------------------------------------ 时间，速度，加速度公式 ------------------------------------------------------------------------------*/ int n,k; HalfKick(); /* 计算半步长内速度 v(t+Dt/2) */ for (n=0; n<nAtom; n++) /* 计算半步长内坐标（距离）s=s+vt */ for (k=0; k<3; k++) r[n][k] = r[n][k] + DeltaT*rv[n][k]; ComputeAccel(); /* 计算加速度, a(t+Dt) */ HalfKick(); /* 计算下半步长速度 v(t+Dt/2) */ } /*----------------------------------------------------------------------------*/ void HalfKick() { /*------------------------------------------------------------------------------ 计算各原子速度rv, v=v0+at. ------------------------------------------------------------------------------*/ int n,k; for (n=0; n<nAtom; n++) for (k=0; k<3; k++) rv[n][k] = rv[n][k] + DeltaTH*ra[n][k]; } /*-----------------