分子动力学模拟的 Matlab 实现.pdf

分子动力学模拟的 Matlab 实现 金剑锋 材料学院 东北大学 氩原子体系分子动力学的 Matlab 实现–主函数 %% 主函数 N = 4*4*4; % 原子总数 r = zeros(N, 3); % 每个原子位置数组 r(x, y, z) v = zeros(N, 3); % 每个原子速度数组 v(x, y, z) a = zeros(N, 3); % 每个原子加速度数组 a(x, y, z) Len = 10; % 模拟盒子的边长 vMax = 0.1; % 速度初始化时的最大值 [r, v] = initialize (r, v, N, Len, vMax); % 初始化位置r 和速度v 0 0 dt = 0.01; % 时间步长dt fileID = fopen(temperture.out,w); 氩原子体系分子动力学的 Matlab 实现–主函数 for time = 1 : 1 : 3000 % 总计算步长 [r, v, a] = velocityVerlet(r, v, a, N, dt); % Verlet 算法 Temp(time) = instantaneousTemperature (v, N); % 计算即时温度 fprintf (Timestep %d Temp %10.2f \n, time, Temp(time)); fprintf (fileID, Timestep %d Temp %10.2f\n, time, Temp(time)); end save(‘data.mat’,‘Temp’); %.mat文件为 matlab标准存储数据格式 fclose(fileID); % 主函数结束 Verlet 子函数的定义 function [r, v, a] = velocityVerlet (r, v, a, N,dt) [a] = computeAccelerations (r, a, N); %计算 t 时刻的加速度 for i =1 : 1 : N for k = 1 : 1 : 3 % r(t+dt) = r(t) + v(t)*dt + 1/2*a(t)*(dt^2) r(i, k) = r(i, k) + v(i, k)*dt + 0.5*a(i, k) *dt*dt ; % v(t+dt) = v(t) + 1/2*a(t)*dt v(i, k) = v(i, k) + 0.5*a(i, k)*dt end end ... end % 子程序结束 计算加速度子函数的定义 function [a] = computeAccelerations(r, a, N) a = zeros(N, 3); % 初始化加速度 a 数组 mass=1.0; eps=1.0; sigma=1.0; %ε eps,σ sigma in LJ potential for i = 1 : 1 : N   14  8  24ε σ σ for j = 1 : 1 : N f(r ) = 2  −     r