24小时热门版块排行榜    

查看: 4856  |  回复: 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的回帖
相关版块跳转 我要订阅楼主 安安安子 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +3 uKOSCocCb6Fm 2026-11-12 6/300 2026-11-15 02:59 by xSrQNSSjcgoX
[论文投稿] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +4 uKOSCocCb6Fm 2026-11-12 10/500 2026-11-15 02:15 by xSrQNSSjcgoX
[文学芳草园] 厉娜 +4 张丹枫1991 2026-11-09 6/300 2026-11-15 00:48 by Equinoxhua
[硕博家园] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +4 G2rhSvqvhNzl 2026-11-10 6/300 2026-11-15 00:27 by Equinoxhua
[考博] 27博士 +6 锦鲤送福至 2026-11-10 6/300 2026-11-14 21:23 by xupeii
[论文投稿] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 XJXMuRjD7cK7 2026-11-13 3/150 2026-11-14 13:02 by KXypaGgBBlw7
[硕博家园] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +6 zSbvVX9JFl9J 2026-11-09 9/450 2026-11-13 11:13 by Equinoxhua
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +6 FChKX8meu1ds 2026-11-09 10/500 2026-11-13 11:08 by Equinoxhua
[硕博家园] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +4 Z7qfCsdIiD1T 2026-11-10 6/300 2026-11-12 11:22 by HajghKOMynmz
[考博] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +4 czAiolzR0IJC 2026-11-10 7/350 2026-11-12 11:11 by HajghKOMynmz
[教师之家] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +4 czAiolzR0IJC 2026-11-10 7/350 2026-11-12 11:00 by HajghKOMynmz
[考研] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +5 C6CRn66RYYXS 2026-11-10 7/350 2026-11-12 10:16 by HajghKOMynmz
[公派出国] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +7 FChKX8meu1ds 2026-11-09 9/450 2026-11-12 10:03 by HajghKOMynmz
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +6 FChKX8meu1ds 2026-11-09 11/550 2026-11-12 09:57 by HajghKOMynmz
[考研] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +6 FChKX8meu1ds 2026-11-09 8/400 2026-11-12 09:41 by HajghKOMynmz
[考研] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 cdqGD6ftHmbl 2026-11-11 3/150 2026-11-12 07:53 by HajghKOMynmz
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 C6CRn66RYYXS 2026-11-10 5/250 2026-11-12 06:03 by 9fOTqj074tIg
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +4 cqv1zKnzd7hr 2026-11-09 8/400 2026-11-12 05:08 by 9fOTqj074tIg
[考博] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +6 zSbvVX9JFl9J 2026-11-09 9/450 2026-11-12 04:57 by 9fOTqj074tIg
[教师之家] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +6 cqv1zKnzd7hr 2026-11-09 7/350 2026-11-12 00:44 by 9fOTqj074tIg
信息提示
请填处理意见