金属铝分子动力学模拟

物理计算与设计报告书
院（系）名称：
学生姓名：
专业名称：
班级：
时间：
金属铝分子动力学模拟
摘要：分子动力学模拟，是指对于原子核和电子所构成的多体系统，用计算机模拟
原子核的运动过程，并从而计算系统的结构和性质，其中每一原子核被视为在全部其它
原子核和电子所提供的经验势场作用下按牛顿定律运动。

我们用c语言编写程序，VMD
动画演示得到原子在拉伸过程中的变化。

在控制温度不变的情况下，得到了金属铝分子
的动力学模拟过程。

通过不断拉伸，趋衡铝分子，计算其势能，力,速度，观察每次拉
伸过程中以及拉伸后铝原子的排列，得到金属铝的运动细节，从而更加利于我们了解铝
的性质。

结论：原子两端的拉力与原子势能的变化曲线基本一致。

原子间断层以滑层方式断
裂。

关键词：铝分子，分子动力学，c语言，势能
1 引言
人们很早就知道材料的力学性能随尺度发生变化尺度减小, 材料中缺陷存在的几率降低, 材料的强度提高同时尺度的变化可能导致材料内在变形竞争机制的改变, 例如多晶材料晶粒粒径在微米级以上时, 强度主要受位错强化机制控制, 而粒径进入纳米级后, 材料的变形主要来源于晶界滑移等机制原子尺度下, 微观效应占主导地位, 材料的理化、力学性能表现出与宏观不同、甚至相反的特性。

Brenner发现金属单晶晶须拉伸强度与晶须直径呈反比,Fleck在微米级细铜丝的扭转试验中观察到尺寸效应纳米电机系统(NEMS)的出现同迫切要求了解纳米尺度下材料的力学行为, 当前从实验上较难获得详细的信息, 而分子动力学模拟可以提供相关细节.
分子动力学通过直接模拟原子的运动过程, 使我们能够详细了解模拟对象的演化发展历史分子动力学模拟的一个关键在于原子势函数的选取原子势早期一般采用简单的对势, 但对势无法正确描述弹性常数, 其结果不理想世纪年代提出的镶嵌原子法、有效介质理论更客观地反映了原子间多体作用的本质, 可得到较合理的结果.认为体系总能量为
2
/102
0]12exp[]1exp[∑
∑∑∑
⎥⎥⎦
⎤
⎢⎢⎣⎡
⎪⎪⎭⎫ ⎝⎛---⎪⎪⎭
⎫ ⎝⎛--=
≠≠i
i
j ij i
i
j ij tol r r q r r p A E ξ
有效地对分子动力学的数据进行后处理也是一个重要的研究方向, 因为对数据进行后处理, 获取有用的信息也是一个很烦琐的工作。

现有的一些分子动力学软件功能不够强大, 只是偏重于某一行业, 通用性不高, 还有一些软件为自由软件, 可维护性不强。

还没有出现集建模、求解、后处理于一体的分子动力学软件。

因此, 针对特定的问题进行自行编程显得尤为重要。

通过分子动力学方法,模拟了铝分子拉伸实验中的形变过程. 研究了晶体取向裂纹的形变特点和断裂机理,观察到各种形变现象,如位错形核和发射,位错运动,堆垛层错或孪晶的形成,纳米空洞的形成与连接等. 计算结果表明,裂纹扩展是塑性过程和弹性过程相结合的过程,其中塑性过程表现为由裂尖发射的位错导致的原子切变行为,而弹性过程的发生则是由无位错区中的原子断键所导致. 本研究采用VC++自行编程, VMD 动画演示得到原子在拉伸过程中的变化,对所研究的问题进行求解。

2 原理
I) 其基本原理是使用一个含有有限个分子并且有周期性边界条件的立方盒子,从该体系某一设定的位能模型出发,通过计算机模拟求解微元中全部分子的牛顿运动方程,记录它们在各个不同时刻的位置、速度和受力等,然后统计得到体系的各种热力学、结构和性质,也就是由体系粒子的微观性质求算其宏观性质。

进行分子动力学模拟的首要问题是要得到准确的原子之间的相互作用势函数。

知道原子间正确的相互作用势,从而必须知道相应的电子基态。

计算中根据以下基本假设:
（1）所有粒子的运动都遵循经典牛顿力学规律。

(2) 粒子之间的相互作用满足叠加原理。

显然这两条忽略了量子效应和多体作用, 与真实物理系统存在一定差别, 仍然属于近似计算。

II)VMD 作图:VMD 是一个强大的原子作图及动画演示软件,在运用C 语言知识对上面求解计算后将得到一个zuobiao.xyz 文件，将此文件直接拖进VMD 原子作图工具便可以得到原子运动动画。

3 方法
第一步: 即模型的设定，也就是势函数的选取。

势函数的研究和物理系统上对物质的描述研究息息相关。

最早是硬球势，即小于临界值时无穷大，大于等于临界值时为零。

常用的是LJ 势函数，还有EAM 势函数，不同的物质状态描述用不同的势函数。

第二步:给定初始条件，也就是给定原子的空间位置和速度。

运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。

如：verlet 算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。

第三步：利用公式：2)()()()(2t dt t dt t dt t i i i i a v r r ++=+计算在第n+1步时所有粒子所处的空间位置。

计算在第n+1步时所有粒子的速度
[]2)()()()(t dt t dt t dt t i i i i a a v v +++=+，动能和速度标度因子：
2
/1212
1)()22(,)(2
1⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡
-==
∑∑
++i n i i
n i k v m kT N v m E β 计算每个时间步系统的总势能U ，考察系统
达到平衡所需要的时间。

第四步：计算将速度乘以标度因子的值，并让该值作为下一次计算时，第n+1步粒子的速度。

第五步：返回第三步，开始第n+2步的模拟计算。

重复几千次，一直到系统达到平衡。

第六步：对系统进行位移加载，每次加载后使系统平衡1000步，让系统经历一个准静态加载的过程。

每次加载并让系统平衡后，计算宏观量。

4 结果
1）程序运行部分结果用VMD 作图的结果：
图（1）
图（2）
图（3）2）程序运行结果，势能用matlab作图的结果：I).利用matlab演示力随应变拉伸的变化：
II).利用matlab演示势能随时间的变化：
图（5）
5 讨论
图（1）图（2）图（3）是铝分子拉伸过程的一组图，分别表示拉伸前，中，后期。

可
以发现在拉伸前铝分子的前两排和后两排的原子保持固定不变，分别向上下两方向拉伸，拉伸一次，趋衡一次并且得到每次原子在平衡状态下的坐标和势能。

由图（3）又
可知，拉伸次数过多导致铝分子断裂。

由图（4）图（5）可发现力的变化曲线与势能的
变化趋于一致，说明势能与力呈一定的相关性。

参考文献
[1]罗熙淳,梁迎春,董申. 单晶铝纳米切削过程分子动力学模拟技术研究.中国机械工程, 2000；11(8):860-8621.
[2]曹莉霞，王崇愚.α2Fe 裂纹的分子动力学研究.钢铁研究总院功能材料所,北京100081
[3]张媛媛, 张文飞, 王鹏。

不同温度下纳米单晶铜杆拉伸的分子动力学模拟燕山大学建筑工程与力学学院, 河北秦皇岛066004
附录：
采用VC++编程:
#include <stdio.h>
#include <math.h>
#include<stdlib.h>
double fen;
double rij;
double A=0.1221,ebsl=1.316,r0=2.86378,wendu=10;
int N=300,m,s,t,n,i,j;
double sm=0.0000861738569;//k
double q=2.516,heli,xheli,eksj;
double li[300][2];
double sheli1,sheli2,p=8.612;
double a[300][2],v[300][2],lii[300][2],shineng[1000],shineng3; double dt=0.5,mal=27,dn,bdyz,kk;
void lih ();//力的函数的声明
void lih ()//计算力
{
for(i=0;i<N;i++)
{
shineng[i]=0;
li[i][0]=0.0;
li[i][1]=0.0;
}
for(i=0;i<N;i++)
{
fen=0.0;
for(j=0;j<N;j++)
{
if(j!=i)
{
rij=pow((a[i][0]-a[j][0]),2)+pow((a[i][1]-a[j][1]),2); rij=sqrt(rij);
fen+=exp(-2*q*(rij/r0-1));
}
}
xheli=sqrt(fen);
if(j!=i)
{
shineng[i]=-ebsl*xheli;
shineng[i]+=A*exp(-p*(rij/r0-1));
}
for(j=0;j<N;j++)
{
if(j!=i)
{
rij=pow((a[i][0]-a[j][0]),2)+pow((a[i][1]-a[j][1]),2);
rij=sqrt(rij);
sheli1=A*exp((-p)*(rij/r0-1))*(p/r0)/rij;
sheli2=ebsl*exp((-2)*q*(rij/r0-1))*(2*q)/r0/rij/xheli;
li[i][0]+=(sheli1-sheli2)*(a[i][0]-a[j][0]);
li[j][0]+=-(sheli1-sheli2)*(a[i][0]-a[j][0]);
li[i][1]+=(sheli1-sheli2)*(a[i][1]-a[j][1]);
li[j][1]+=-(sheli1-sheli2)*(a[i][1]-a[j][1]);
}
}
}
}
void main()//主函数
{
double a1=sqrt(2)*4.050/2,a2=sqrt(3/2)*4.050/2;
int i,j,k=0;
FILE*fp;
for (i=0;i<20;i++)
{
for (j=0;j<15;j++)
{
if (i%2==0)
{
a[k][0]=j*a1;
a[k][1]=a2*i;
}
else
{
a[k][0]=j*a1+a1/2;
a[k][1]=a2*i;
}
k++;
}
}
fp=fopen("zuobiao.dat","w");
fprintf(fp,"\n\n");
for (i=0;i<k;i++)
{
fprintf(fp,"Al %.18f %.18f %.12f\n",a[i][0],a[i][1],0.0); }
fclose(fp);
//初始速度
for(i=0;i<N;i++)
{
v[i][0]=0.0;
v[i][1]=0.0;
}
lih (); //调用函数
//输出势能
fp=fopen("shineng.dat","w");
fprintf(fp,"\n\n");
fprintf(fp,"AL 势能=%.16f\n",shineng[N]);
fclose(fp);
//趋衡过程
for (s=0;s<100;s++)//趋衡100次
{
for(i=0;i<N;i++)
{
a[i][0]=a[i][0]+dt*v[i][0]+0.5*dt*dt*li[i][0]/mal; a[i][1]=a[i][1]+dt*v[i][1]+0.5*dt*dt*li[i][1]/mal; lii[i][0]=li[i][0];
lii[i][1]=li[i][1];
li[i][0]=0.0;
li[i][1]=0.0;
}
lih(); //重新计算力
//算速度
for(i=0;i<k;i++)
{
v[i][0]+=0.5*dt*(li[i][0]+lii[i][0])/mal;
v[i][1]+=0.5*dt*(li[i][1]+lii[i][1])/mal;
}
dn=0.0;
for(n=0;n<N;n++)
{
dn+=0.5*mal*(v[n][0]*v[n][0]+v[n][1]*v[n][1]);
}
//标度因子
bdyz=sqrt(0.5*(2*n-3)*wendu*sm/dn);
//对每个原子的速度进行标度
for(n=0;n<N;n++)
{
v[n][0]=v[n][0]*bdyz;
v[n][1]=v[n][1]*bdyz;
}
}
fp=fopen("zuobiao.dat","a");
fprintf(fp,"\n\n");
for(i=0;i<k;i++)
{
fprintf(fp,"AL %.16f %.16f %.16f\n",a[i][0],a[i][1],0.0); }
fclose(fp);
//拉伸阶段
for(t=0;t<500;t++)//拉伸500次
{
for(i=N-30;i<N;i++)
{
a[i][1]=a[i][1]+0.0005*15*a2;
}
for(i=0;i<30;i++)
{
a[i][1]=a[i][1]-0.0005*15*a2;
}
//趋衡过程
for (s=0;s<50;s++)//趋衡50次
{
for(i=30;i<N-30;i++)
{
a[i][0]=a[i][0]+dt*v[i][0]+0.5*dt*dt*li[i][0]/mal;
a[i][1]=a[i][1]+dt*v[i][1]+0.5*dt*dt*li[i][1]/mal;
lii[i][0]=li[i][0];
lii[i][1]=li[i][1];
li[i][0]=0.0;
li[i][1]=0.0;
}
//重新计算力
for(i=0;i<N;i++)
{
shineng3=0.0;
fen=0.0;
for(j=0;j<N;j++)
{
if(j!=i)
{
rij=pow((a[i][0]-a[j][0]),2)+pow((a[i][1]-a[j][1]),2); rij=sqrt(rij);
fen+=exp(-2*q*(rij/r0-1));
}
}
xheli=sqrt(fen);
if(j!=i)
{
shineng3=-ebsl*xheli;
shineng3+=A*exp(-p*(rij/r0-1));
}
for(j=0;j<N;j++)
{
if(j!=i)
{
rij=pow((a[i][0]-a[j][0]),2)+pow((a[i][1]-a[j][1]),2);
rij=sqrt(rij);
sheli1=A*exp((-p)*(rij/r0-1))*(p/r0)/rij;
sheli2=ebsl*exp((-2)*q*(rij/r0-1))*(2*q)/r0/rij/xheli;
li[i][0]+=(sheli1-sheli2)*(a[i][0]-a[j][0]);
li[j][0]+=-(sheli1-sheli2)*(a[i][0]-a[j][0]);
li[i][1]+=(sheli1-sheli2)*(a[i][1]-a[j][1]);
li[j][1]+=-(sheli1-sheli2)*(a[i][1]-a[j][1]);
}
}
}
shineng[t]=shineng3;
//算速度
for(i=0;i<k;i++)
{
v[i][0]+=0.5*dt*(li[i][0]+lii[i][0])/mal;
v[i][1]+=0.5*dt*(li[i][1]+lii[i][1])/mal;
}
dn=0.0;
for(n=0;n<N;n++)
{
dn+=0.5*mal*(v[n][0]*v[n][0]+v[n][1]*v[n][1]);
}
//标度因子
bdyz=sqrt(0.5*(2*n-3)*wendu*sm/dn);
//对每个原子的速度进行标度：
for(n=0;n<N;n++)
{
v[n][0]=v[n][0]*bdyz;
v[n][1]=v[n][1]*bdyz;
}
}
//输出坐标
printf ("第%d次",t+1);
fp=fopen("zuobiao.dat","a");
fprintf(fp,"\n\n");
for(i=0;i<k;i++)
{
fprintf(fp,"AL %.16f %.16f %.16f\n",a[i][0],a[i][1],0.0); }
fclose(fp);
//输出势能
fp=fopen("shineng.dat","a");
fprintf(fp,"\n\n");
fprintf(fp,"AL 势能=%.16f\n",shineng[t]);
fclose(fp);
}//拉伸的括号
}//main函数的。