两原子的分子动力学模拟-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,:);%将质心速度设置为0

d=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);%总能

for i=1:n-1

r1(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);

matlab 分子动力学,两体的分子动力学模型-MATLAB源程序相关推荐

  1. matlab中两个符号矩阵相加,MATLAB矩阵分析及符号运算.ppt

    MATLAB矩阵分析及符号运算 第三讲 MATLAB的符号运算 -- matlab 不仅具有数值运算功能,还开发了在matlab环境下实现符号计算的工具包Symbolic Math Toolbox 符 ...

  2. matlab stract结构_6. matlab入门——结构体、元胞数组、字符串

    1.结构体 (1)使用赋值方法创建结构体 %% 使用赋值方法创建结构体 person(1).name = '张三'; person(1).weight = 66; person(1).length = ...

  3. 【 MATLAB 】两个序列的卷积和运算的MATLAB实现(2)

    已知下面两个序列: 求这两个序列的卷积. 求卷积的函数是conv,但是使用这个函数有个问题,就是下标问题,也就是求卷积之后的元素值的位置.因此,我们必须要定一个起始点和一个结束点. 方法: 是两个有限 ...

  4. matlab if语句多个执行举例,初学Matlab,有两个语句,if语句和switch语句,有两个例子哪位大神能帮我讲讲...

    问题描述: 初学Matlab,有两个语句,if语句和switch语句,有两个例子哪位大神能帮我讲讲 if logical_expression statements elseif logical_ex ...

  5. matlab根据结构体数组,用邻接矩阵和序遍历创建树形结构:

    matlab根据结构体数组,用邻接矩阵和先序遍历创建树形结构: https://blog.csdn.net/C_Redrock/article/details/84980241

  6. matlab 判断两个矩阵有元素相等_如何使用MATLAB对Excel中的多参数进行计算?

    THE START MATLAB和Excel这两者之间有着什么样的关系呢?今天我把之前学习以及用到的关于用MATLAB读写Excel数据,并进行计算处理的经验分享给需要的小伙伴.参加过数学建模的这个应 ...

  7. matlab struct 结构体

    matlab 的结构体第一次让 matlab 中的变量有了可以通过 .访问的成员变量,有了类的含义,甚至是面向对象的意味. 1. 结构体的赋值 结构体的赋值,这里不建议用下面这种形式进行统一赋值, s ...

  8. matlab 两个数中取小,matlab中取两个数中的较小值

    在EXCEL表中,如何取一组数据中的两个最大数和两个最小数?用什么函数? =large(a1:a15,1)第一大=large(a1:a15,2)第二大=small(a1:a15,1)第一小=small ...

  9. centos 7安装matlab的两种方法(桌面安装和命令行安装)

    matlab安装说明 安装之前一直以为命令行安装(静默安装)完就是命令行界面,安装成功后才发现还是有桌面的,还跟桌面安装的一模一样.所以,个人建议对linux不太熟悉的还是用桌面版安装,虽然会有点卡顿 ...

  10. matlab 求曲面体积,matlab求两曲面之间的体积

    MATLAB求曲面相交所成空间曲线的图形 放在你程序后也可,单独运行也行:t=-0.1:0.1:2*pi;x=2*cos(t);%交线参数方程z=2*sin(t);y1=sqrt(5)*ones(si ...

最新文章

  1. 分享篇--esp32直连天猫精灵
  2. C# GDI+ 画坐标(x,y)
  3. C语言基础:C语言指针(6) - 指针和字符串
  4. Scala的控制结构
  5. UnimplementedError: Fused conv implementation does not support grouped convolutions for now
  6. 药店管理系统/APP/小程序/网站
  7. win7 升级到 win10
  8. Installer User Interface Mode Not Supported解决方法
  9. MySQL数据分析实战-朱元禄-专题视频课程
  10. nginx中的timeout超时设置,请求超时、响应等待超时等
  11. serviceWorker 服务器与浏览器之间的代理
  12. NVIDIA Jetson官网资料整理
  13. 基于Docker搭建DzzOffice与OnlyOffice线上协同办公服务器
  14. maven基础:mvn命令常用参数整理;如:-am构建指定模块,同时构建指定模块依赖的其他模块
  15. MNE-Python读取MATLAB保存的.mat文件
  16. springcloud 项目maven依赖:Failure to find org.springframework.cloud:spring-cloud-dependencies
  17. openvswitch 2.3.1 配置详解
  18. java:IO流(缓冲流、对象流、控制台IO、转换流、java.io.File 类 )
  19. c语言程序综合实习学生成绩,C语言程序设计综合实习报告-资源下载人人文库网...
  20. 【仿真】Carla之收集数据快速教程 (附完整代码)

热门文章

  1. 王之泰201771010131《面向对象程序设计(java)》第四周学习总结
  2. 计算机用户界面的设计,计算机软件用户界面设计的基本原则
  3. C语言中的指针加减偏移量
  4. 动手学深度学习(二十七)——微调(fine turning)
  5. 活化脂修饰NOTA,NOTA-NHS ester,CAS:1338231-09-6
  6. 【框架-MFC】MFC 显示和隐藏 星号密码 以及如何预防被查看
  7. [经典]PK:星际争霸 vs 魔兽争霸3 vs 红警
  8. 在word中公式太长,用公式编辑器怎样设置才能自动换行?
  9. 未转变者服务器关雨指令,Unturned未转变者3.21版本物品ID代码汇总
  10. java 数学公式编辑器_妈妈再也不用担心我的公式写不出来了:一款公式输入神器实测...