24小时热门版块排行榜    

查看: 4312  |  回复: 8
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

安安安子

新虫 (初入文坛)

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

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

» 猜你喜欢

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

已阅   回复此楼   关注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的回帖
查看全部 9 个回答

chuntao118

新虫 (著名写手)

【答案】应助回帖

感谢参与,应助指数 +1
可以参考这本书《Understanding Molecular Simulation From Algorithms to Applications》
2楼2016-03-22 14:32:44
已阅   回复此楼   关注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的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 305分求调剂(食品工程) +4 Sxy112 2026-03-21 6/300 2026-03-22 15:26 by 无懈可击111
[考研] 一志愿中南大学化学学硕0703总分337求调剂 +3 niko- 2026-03-22 3/150 2026-03-22 15:15 by 杨杨杨紫
[考研] 260求调剂 +3 朱芷琳 2026-03-20 4/200 2026-03-22 15:12 by 朱芷琳
[考研] 0703化学调剂 +4 妮妮ninicgb 2026-03-21 4/200 2026-03-21 18:39 by 学员8dgXkO
[考研] 工科0856求调剂 +3 沐析汀汀 2026-03-21 3/150 2026-03-21 18:30 by 学员8dgXkO
[考研] 277材料科学与工程080500求调剂 +6 自由煎饼果子 2026-03-16 6/300 2026-03-21 17:21 by 学员8dgXkO
[考研] 268求调剂 +9 简单点0 2026-03-17 9/450 2026-03-21 15:37 by lature00
[考研] 324分 085600材料化工求调剂 +4 llllkkkhh 2026-03-18 4/200 2026-03-21 01:24 by JourneyLucky
[考研] 271材料工程求调剂 +8 .6lL 2026-03-18 8/400 2026-03-21 00:58 by JourneyLucky
[考研] 294求调剂材料与化工专硕 +15 陌の森林 2026-03-18 15/750 2026-03-20 23:28 by JourneyLucky
[考研] 一志愿吉林大学材料学硕321求调剂 +11 Ymlll 2026-03-18 15/750 2026-03-20 19:40 by 丁丁*
[考研] 086500 325 求调剂 +3 领带小熊 2026-03-19 3/150 2026-03-20 18:38 by 尽舜尧1
[考研] 320求调剂0856 +3 不想起名字112 2026-03-19 3/150 2026-03-19 22:53 by 学员8dgXkO
[考研] 材料与化工求调剂 +7 为学666 2026-03-16 7/350 2026-03-19 14:48 by 尽舜尧1
[考研] 085601专硕,总分342求调剂,地区不限 +5 share_joy 2026-03-16 5/250 2026-03-18 14:48 by haxia
[考研] 308求调剂 +4 是Lupa啊 2026-03-16 4/200 2026-03-17 17:12 by ruiyingmiao
[考研] 085601求调剂 +4 Du.11 2026-03-16 4/200 2026-03-17 17:08 by ruiyingmiao
[考研] 290求调剂 +3 p asserby. 2026-03-15 4/200 2026-03-17 16:35 by wangkm
[考研] 有没有道铁/土木的想调剂南林,给自己招师弟中~ +3 TqlXswl 2026-03-16 7/350 2026-03-17 15:23 by TqlXswl
[考研] 070303 总分349求调剂 +3 LJY9966 2026-03-15 5/250 2026-03-16 14:24 by xwxstudy
信息提示
请填处理意见