24小时热门版块排行榜    

查看: 4165  |  回复: 8

安安安子

新虫 (初入文坛)

[求助] 从LAMMPS得到的原子速度,如何计算速度自相关函数(VAF)?求助 已有3人参与

比如说我的系统里有4000个原子,在某个状态下跑了100步,每一步的时间间隔是确定的(2fs),每一步时每个原子的三维速度都会输出。我想计算这个系统的速度自相关函数VAF,计算公式是<v(0)·v(t)> 。我不太明白这个公式的含义是计算什么,用FORTRAN或MATLAB应该如何编程实现呢?紧急求助,不胜感激。
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

chuntao118

新虫 (著名写手)

【答案】应助回帖

感谢参与,应助指数 +1
可以参考这本书《Understanding Molecular Simulation From Algorithms to Applications》
2楼2016-03-22 14:32:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

challenger9

银虫 (初入文坛)

引用回帖:
2楼: Originally posted by chuntao118 at 2016-03-22 14:32:44
可以参考这本书《Understanding Molecular Simulation From Algorithms to Applications》

这个问题你应该早就破解了吧!!我现在在做苯的自相关,求助,求分享
3楼2016-06-24 10:06:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

brucefan

专家顾问 (著名写手)

【答案】应助回帖

% a matlab script for calculating VACF from velocity data
clear;
load v.txt;
% assume your velocity data are in the above file and assume the format is (N is the number of atoms, M is the number of time points):
% vx_1 vy_1 vz_1      % time point 1
% vx_2 vy_2 vz_2      % time point 1
% ...                            % time point 1
% vx_N, vy_N, vz_N   % time point 1
% vx_1 vy_1 vz_1      % time point 2
% vx_2 vy_2 vz_2      % time point 2
% ...                            % time point 2
% vx_N, vy_N, vz_N  % time point 2
% ...                           
% vx_1 vy_1 vz_1     % time point M
% vx_2 vy_2 vz_2     % time point M
% ...                           % time point M
% vx_N, vy_N, vz_N  % time point M

N = xxx; % number of atoms in your system
M = length(v)/N; % number of time points for your velocity data
dt = xxx; % the time interval between two set of velocities (in some unit)
Nt = xxx; % maximum length of the correlation your want (usually Nt = M/10 is a good choice)
time = dt*(0:Nt-1);
M = M-Nt; % you have to waste a small portion of data

vacf=zeros(Nt,1);
for nt=0:Nt-1
   for m=1:M
       vacf(nt+1,=sum(sum(v_all((m-1)*N+1:m*N,.*v_all((m+nt-1)*N+1m+nt)*N,));
   end
end
%vacf=vacf/M;  % you can also normalize it by using [vacf=vacf(1);] if you want

% now you can plot the result:
close all;
figure;
plot(time, vacf,'o-');
xlabel('time (some unit)');
ylabel('VACF (some unit)');
4楼2016-06-27 18:45:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

brucefan

专家顾问 (著名写手)

【答案】应助回帖

sorry, v_all should be v in some line of the code.
5楼2016-06-27 18:52:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cuso4852cuoh

新虫 (小有名气)

【答案】应助回帖

求助自相关函数和声子态密度啊?求楼主
6楼2018-11-26 19:15:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

looyady

金虫 (小有名气)

引用回帖:
6楼: Originally posted by cuso4852cuoh at 2018-11-26 19:15:36
求助自相关函数和声子态密度啊?求楼主

同求声子态密度
7楼2018-12-14 15:21:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

流落的小猫

新虫 (初入文坛)

这个我做过,主要是从lammps中导出数据,然后自己编写程序来做,我通过计算得到的数据曲线还是比较满意的,仅供楼主参考,这是一种方法。
8楼2018-12-19 16:27:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

三岁薇恩

新虫 (正式写手)

引用回帖:
6楼: Originally posted by cuso4852cuoh at 2018-11-26 19:15:36
求助自相关函数和声子态密度啊?求楼主

求到了吗 兄弟
9楼2019-04-10 15:39:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 安安安子 的主题更新
信息提示
请填处理意见