两体的分子动力学模型-MATLAB源程序-

1、两原子的分子动力学模拟-MATLAB程序%假设两个原子的质量相同,都设为1%两原子之间的相互作用势为Lennard-Jones势clear;clc;n=1000;%模拟步数,可修改r1=zeros(n,3;%原子1的坐标r2=zeros(n,3;%原子2的坐标v1=zeros(n,3;%原子1的速度v2=zeros(n,3;%原子2的速度f1=zeros(n,3;%原子1受到的力f2=zeros(n,3;%原子2受到的力kinetic=zeros(n,1;%动能potential=zeros(n,1;%势能total=zeros(n,1;%总能量r1(1,1=0.2*rand(+0.3;r1(

2、1,2=-0.2*rand(-0.1;r1(1,3=0.1*rand(-0.3;%随机设置原子1的初始位置r2(1,:=-r1(1,:;%将系统的质心设置在原点v1(1,:=rand(1,3;%速度服从标准正态分布v2(1,:=-v1(1,:;%将质心速度设置为0d=sqrt(r1(1,1-r2(1,12+(r1(1,2-r2(1,22+(r1(1,3-r2(1,32;%初始时刻两原子间的距离k=48*(d(-14-0.5*(d(-8;%力关于坐标的比例系数f1(1,:=(r1(1,:-r2(1,:*k;%初始时刻的受力f2(1,:=-f1(1,:;%牛顿第三定律h=0.01;%模拟步长,步长

3、过长会出现机械能不守恒的情况kinetic(1,1=v1(1,12+v1(1,22+v1(1,32;%初始动能,两原子具有相同的动能potential(1,1=4*(d(-12-d(-6;%初始势能total(1,1=kinetic(1,1+potential(1,1;%总能for i=1:n-1r1(i+1,:=r1(i,:+v1(i,:*h+0.5*f1(i,:*h2;r2(i+1,:=-r1(i+1,:;d=sqrt(r1(i+1,1-r2(i+1,12+(r1(i+1,2-r2(i+1,22+(r1(i+1,3-r2(i+1,32;k=48*(d(-14-0.5*(d(-8;f1(i+

4、1,:=(r1(i+1,:-r2(i+1,:*k;f2(i+1,:=-f1(i+1,:;v1(i+1,:=v1(i,:+0.5*h*(f1(i+1,:+f1(i,:;v2(i+1,:=-v1(i+1,:;kinetic(i+1,1=(v1(i+1,12+v1(i+1,22+v1(i+1,32;%两原子的动能相同potential(i+1,1=4*(d(-12-d(-6;total(i+1,1=kinetic(i+1+potential(i+1;endt=h:h:n*h;%模拟时间subplot(121,plot(t,kinetic,'g',t,potential,'r',t,total,'b'legend('动能','势能','总能'xlabel('时间'ylabel('能量'title('能量-时间变化曲线'grid on;subplot(122,plot3(r1(:,1,r1(:,2,r1(:,3,'b',r2(:,1,r2(:,2,r2(:,3,'g' legend('原子1','原子2'grid on;xlabel('x'ylabel(