24小时热门版块排行榜    

Znn3bq.jpeg
查看: 4410  |  回复: 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的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 26药学专硕105500求调剂 +6 喽哈加油 2026-04-13 6/300 2026-04-14 16:40 by zhouxiaoyu
[考研] 恳请有学校收留 +3 柯淮然 2026-04-12 3/150 2026-04-14 16:25 by 逆水乘风
[考研] 071000生物学调剂求助 +18 zzzzwww 2026-04-09 21/1050 2026-04-14 15:39 by zs92450
[考研] 本科211,报考085601-310分 +16 ararak 2026-04-13 16/800 2026-04-14 14:55 by Delta2012
[考研] 求调剂 +16 雪逢冬 2026-04-10 16/800 2026-04-14 14:27 by 逆水乘风
[考研] 一志愿双非085400电子信息344 求调剂,对材料和化学方向也感兴趣 +12 无情的小羊 2026-04-09 13/650 2026-04-13 14:17 by 张zhihao
[考研] 求调剂288 +7 ioodiiij 2026-04-10 9/450 2026-04-13 08:33 by Hayaay
[考研] 346分,工科0854求调剂,专硕 +6 moser233 2026-04-12 7/350 2026-04-12 22:11 by fqwang
[考研] 本科南方医科大学 一志愿985 药学学硕284分 求调剂 +5 弱水听文 2026-04-09 5/250 2026-04-12 13:16 by shengxi123
[考研] 电气专硕320求调剂 +6 小麻子111 2026-04-10 6/300 2026-04-12 10:54 by lemon6009
[考研] 296求调剂 +14 汪!?! 2026-04-10 16/800 2026-04-12 10:48 by zhouyuwinner
[考研] 本人女孩 +7 吼吼, 2026-04-10 9/450 2026-04-11 14:45 by ACS Nano——
[考研] 085600材料与化工329分求调剂 +16 叶zilin 2026-04-10 16/800 2026-04-11 11:04 by may_新宇
[考研] 求调剂 +6 archer.. 2026-04-09 8/400 2026-04-11 10:55 by zhq0425
[考研] 工科273调剂 +6 X1999 2026-04-09 7/350 2026-04-11 10:23 by zhq0425
[考研] 中药学调剂 初试324 +4 洋甘菊、 2026-04-10 6/300 2026-04-11 09:41 by gong120082
[考研] 346,工科0854求调剂,专硕 +7 moser233 2026-04-10 8/400 2026-04-11 08:52 by 猪会飞
[考研] 263能源动力专硕求调剂 +4 加大号饭盒袋 2026-04-10 4/200 2026-04-10 20:52 by gong120082
[考研] 一志愿华东师范生物学326分,求调剂 +8 刘墨墨 2026-04-09 8/400 2026-04-10 12:00 by pengliang8036
[考研] 070300化学 求调剂 +13 73372112 2026-04-08 13/650 2026-04-09 20:22 by maddjdld
信息提示
请填处理意见