搜索
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(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,1))^2+(r1(1,2)-r2(1,2))^2+(r1(1,3)-r2(1,3))^2);%初始时刻两原子间的距离k=48*(d^(-14)-0.5*(d^(-8)));%力关于坐标的比例系数f1(1,:)=(r1(1,:)-r2(1,:))*k;%初始时刻的受力f2(1,:)=-f1(1,:);%牛顿第三定律h=0.01;%模拟步长，步长过长会出现机械能不守恒的情况kinetic(1,1)=v1(1,1)^2+v1(1,2)^2+v1(1,3)^2;%初始动能,两原子具有相同的动能potential(1,1)=4*(d^(-12)-d^(-6));%初始势能total(1,1)=kinetic(1,1)+potential(1,1);%总能fori=1:n-1r1(i+1,:)=r1(i,:)+v1(i,:)*h+0.5*f1(i,:)*h^2;r2(i+1,:)=-r1(i+1,:);d=sqrt((r1(i+1,1)-r2(i+1,1))^2+(r1(i+1,2)-r2(i+1,2))^2+(r1(i+1,3)-r2(i+1,3))^2);k=48*(d^(-14)-0.5*(d^(-8)));f1(i+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,1)^2+v1(i+1,2)^2+v1(i+1,3)^2);%两原子的动能相同potential(i+1,1)=4*(d^(-12)-d^(-6));total(i+1,1)=kinetic(i+1)+potential(i+1);2endt=h:h:n*h;%模拟时间subplot(121),plot(t,kinetic,'g',t,potential,'r',t,total,'b');legend('动能','势能','总能');xlabel('时间');ylabel('能量');title('能量-时间变化曲线');gridon;subplot(122),plot3(r1(:,1),r1(:,2),r1(:,3),'b',r2(:,1),r2(:,2),r2(:,3),'g');legend('原子1','原子2');gridon;xlabel('x');ylabel('y');zlabel('z');title('两原子轨迹图');运行结果如图所示（结果具有随机性，本文档附上一种好看的图形）：3