1、兩原子的分子動力學(xué)模擬-MATLAB程序%假設(shè)兩個原子的質(zhì)量相同,都設(shè)為1%兩原子之間的相互作用勢為Lennard-Jones勢clear;clc;n=1000;%模擬步數(shù),可修改r1=zeros(n,3;%原子1的坐標(biāo)r2=zeros(n,3;%原子2的坐標(biāo)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;%隨機設(shè)置原子1的初始位置r2(1,:=-r1(1,:;%將系統(tǒng)的質(zhì)心設(shè)置在原點v1(1,:=rand(1,3;%速度服從標(biāo)準(zhǔn)正態(tài)分布v2(1,:=-v1(1,:;%將質(zhì)心速度設(shè)置為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;%力關(guān)于坐標(biāo)的比例系數(shù)f1(1,:=(r1(1,:-r2(1,:*k;%初始時刻的受力f2(1,:=-f1(1,:;%牛頓第三定律h=0.01;%模擬步長,步長

3、過長會出現(xiàn)機械能不守恒的情況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(