兩體的分子動(dòng)力學(xué)模型-MATLAB源程序-

1、兩原子的分子動(dòng)力學(xué)模擬-MATLAB程序%假設(shè)兩個(gè)原子的質(zhì)量相同,都設(shè)為1%兩原子之間的相互作用勢(shì)為L(zhǎng)ennard-Jones勢(shì)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的速度f(wàn)1=zeros(n,3;%原子1受到的力f2=zeros(n,3;%原子2受到的力kinetic=zeros(n,1;%動(dòng)能potential=zeros(n,1;%勢(shì)能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;%隨機(jī)設(shè)置原子1的初始位置r2(1,:=-r1(1,:;%將系統(tǒng)的質(zhì)心設(shè)置在原點(diǎn)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;%初始時(shí)刻兩原子間的距離k=48*(d(-14-0.5*(d(-8;%力關(guān)于坐標(biāo)的比例系數(shù)f1(1,:=(r1(1,:-r2(1,:*k;%初始時(shí)刻的受力f2(1,:=-f1(1,:;%牛頓第三定律h=0.01;%模擬步長(zhǎng),步長(zhǎng)

3、過(guò)長(zhǎng)會(huì)出現(xiàn)機(jī)械能不守恒的情況kinetic(1,1=v1(1,12+v1(1,22+v1(1,32;%初始動(dòng)能,兩原子具有相同的動(dòng)能potential(1,1=4*(d(-12-d(-6;%初始勢(shì)能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;%兩原子的動(dòng)能相同potential(i+1,1=4*(d(-12-d(-6;total(i+1,1=kinetic(i+1+potential(i+1;endt=h:h:n*h;%模擬時(shí)間subplot(121,plot(t,kinetic,'g',t,potential,'r',t,total,'b'legend('動(dòng)能','勢(shì)能','總能'xlabel('時(shí)間'ylabel('能量'title('能量-時(shí)間變化曲線'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(