poster
2019-11-29, 11:00
我如何将这个MATLAB仿真代码更改为euler算法以进行时间积分？我正在matlab中进行分子动力学模拟。下面是Verlet算法的模拟代码。它是关于伦纳德·琼斯在sigma = 1中进行力计算的潜力。温度是300K，玻尔兹曼常数是1。我想将此Verlet算法代码更改为Euler算法以进行时间积分。
clear all; npart = 100; x = zeros(npart, 1); v = zeros(npart, 1); xm = zeros(npart, 1); f = zeros(npart, 1); %tmax = 1 ; tep = 300; sigma= 1; % dt = time step size %dt = tmax/total_time_step; dt = 0.00001; tmax = 1; % time step (integer) total_time_step = uint8(tmax / dt); temp = zeros(total_time_step,1); pot = zeros(total_time_step,1); kin = zeros(total_time_step,1); etot = zeros(total_time_step,1); time = zeros(total_time_step,1); % Cutoff radius rc=2.56*sigma; %Initial particle distance % ( force equilibrium at r0 = sigma*(2)^(1/6) in LJ potential) eps=1.122462048*sigma; % Total size of particles box=(npart)*eps; t=0; tstep=0; sumv=0; sumv2=0; %initialize positions, velocities of particles for i=1: npart x(i)=(i-1)*eps; % position of 0th time step v(i)=(rand()-0.5)*2; sumv=sumv+v(i); sumv2=sumv2+v(i)^2; end sumv=sumv/npart; sumv2=sumv2/npart; fs=sqrt(tep/sumv2); for i=1: npart v(i)=(v(i)-sumv)*fs; % velocity of 0th time step % xm is the position of particles at step -1 (for verlet algorithm) xm(i)=x(i)-v(i)*dt; % velocity of -1th time step (for verlet algorithm) end xmin = x(1)-0.5*eps; xmax = x(npart)+0.5*eps; while (t < tmax) % nth time step tstep = tstep + 1; time(tstep) = t; % calculate interacting forces for all at nth time step f=zeros(npart, 1); en =0; for i=1: npart-1 for j=i+1: npart xr=x(i)-x(j); xr=xr-box*round(xr/box); r2=xr^2; rc2=rc^2; if (r2 > 0 && r2 < rc2) r=abs(xr); ri=sigma/r; ff = 24*(2*ri^12-ri^6)/(r^2); f(i)=f(i)+ff*xr; f(j)=f(j)-ff*xr; en=en+4*ri^6*(ri^6-1); end end end sumv=0; sumv2=0; % time integration to update position(x) (verlet algorithm) for i=1: npart x(i)= x(i) + dt*v(i) + (dt^2)/2*f(i); % position of (n+1)th time step v(i)= v(i) + dt*f(i); % velocity of (n)th time step sumv=sumv+v(i); sumv2=sumv2+v(i)^2; % for periodic boundary condition if x(i)xmax x(i) = x(i) - box; xm(i) = xm(i) - box; end end temp(tstep)=sumv2/(npart); pot(tstep) =en/npart; kin(tstep)=(0.5*sumv2)/npart; etot(tstep)=(en+0.5*sumv2)/npart; t=t+dt; tstep end figure('name', 'result plot'); subplot(121); plot(time, pot, '-b', 'linewidth', 1); xlabel('time'); ylabel('energy'); hold on; plot(time, etot, '-r', 'linewidth', 1); plot(time, kin, '-g', 'linewidth', 1); subplot(122); plot(time,temp, 'linewidth', 1); xlabel('time'); ylabel('temperature');
更多&回答... (https://stackoverflow.com/q/59098169)
clear all; npart = 100; x = zeros(npart, 1); v = zeros(npart, 1); xm = zeros(npart, 1); f = zeros(npart, 1); %tmax = 1 ; tep = 300; sigma= 1; % dt = time step size %dt = tmax/total_time_step; dt = 0.00001; tmax = 1; % time step (integer) total_time_step = uint8(tmax / dt); temp = zeros(total_time_step,1); pot = zeros(total_time_step,1); kin = zeros(total_time_step,1); etot = zeros(total_time_step,1); time = zeros(total_time_step,1); % Cutoff radius rc=2.56*sigma; %Initial particle distance % ( force equilibrium at r0 = sigma*(2)^(1/6) in LJ potential) eps=1.122462048*sigma; % Total size of particles box=(npart)*eps; t=0; tstep=0; sumv=0; sumv2=0; %initialize positions, velocities of particles for i=1: npart x(i)=(i-1)*eps; % position of 0th time step v(i)=(rand()-0.5)*2; sumv=sumv+v(i); sumv2=sumv2+v(i)^2; end sumv=sumv/npart; sumv2=sumv2/npart; fs=sqrt(tep/sumv2); for i=1: npart v(i)=(v(i)-sumv)*fs; % velocity of 0th time step % xm is the position of particles at step -1 (for verlet algorithm) xm(i)=x(i)-v(i)*dt; % velocity of -1th time step (for verlet algorithm) end xmin = x(1)-0.5*eps; xmax = x(npart)+0.5*eps; while (t < tmax) % nth time step tstep = tstep + 1; time(tstep) = t; % calculate interacting forces for all at nth time step f=zeros(npart, 1); en =0; for i=1: npart-1 for j=i+1: npart xr=x(i)-x(j); xr=xr-box*round(xr/box); r2=xr^2; rc2=rc^2; if (r2 > 0 && r2 < rc2) r=abs(xr); ri=sigma/r; ff = 24*(2*ri^12-ri^6)/(r^2); f(i)=f(i)+ff*xr; f(j)=f(j)-ff*xr; en=en+4*ri^6*(ri^6-1); end end end sumv=0; sumv2=0; % time integration to update position(x) (verlet algorithm) for i=1: npart x(i)= x(i) + dt*v(i) + (dt^2)/2*f(i); % position of (n+1)th time step v(i)= v(i) + dt*f(i); % velocity of (n)th time step sumv=sumv+v(i); sumv2=sumv2+v(i)^2; % for periodic boundary condition if x(i)xmax x(i) = x(i) - box; xm(i) = xm(i) - box; end end temp(tstep)=sumv2/(npart); pot(tstep) =en/npart; kin(tstep)=(0.5*sumv2)/npart; etot(tstep)=(en+0.5*sumv2)/npart; t=t+dt; tstep end figure('name', 'result plot'); subplot(121); plot(time, pot, '-b', 'linewidth', 1); xlabel('time'); ylabel('energy'); hold on; plot(time, etot, '-r', 'linewidth', 1); plot(time, kin, '-g', 'linewidth', 1); subplot(122); plot(time,temp, 'linewidth', 1); xlabel('time'); ylabel('temperature');
更多&回答... (https://stackoverflow.com/q/59098169)