24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

安安安子

新虫 (初入文坛)

应助: 0 (幼儿园)
帖子: 2
在线: 32分钟
虫号: 2410151
注册: 2013-04-09

[求助] 从LAMMPS得到的原子速度，如何计算速度自相关函数（VAF）？求助已有3人参与

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

回复此楼

» 猜你喜欢

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

lammps中计算弹性常数的例子是基于什么原理？已经有1人回复
从lammps模拟使用的势函数画出势能曲线已经有0人回复
急求解答：lammps安装问题已经有10人回复
对原子数过多的晶格进行结构优化已经有3人回复
【转】（已搜索无重复）LAMMPS中fix ave spatial空间平均命令中文版已经有26人回复
求助？怎么从剑桥晶体库得到物质的结构数据已经有5人回复
lammps 时间关联函数求助已经有1人回复
LAMMPS温度－步长曲线震荡已经有11人回复
有用lammps模拟过聚合物的相关性能吗？已经有12人回复
python之蟒开启理论计算化学的新时代已经有22人回复
关于towhee 回答zhang_jaj的问题已经有44人回复
怎样用lammps得到无定形Si的结构已经有3人回复
站里关于MS Analysis计算材料的应力应变与Cubic以及Trigonal cell弹性常数的讨论已经有34人回复

1楼 2016-03-22 14:09:17

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

brucefan

专家顾问 (著名写手)

专家经验: +218
模拟EPI: 4
应助: 422 (硕士)
金币: 2294.8
散金: 1720
红花: 88
帖子: 1423
在线: 670.2小时
虫号: 3048965
注册: 2014-03-14
性别: GG
专业: 凝聚态物性 II ：电子结构
管辖: 分子模拟

【答案】应助回帖

% 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+1

m+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)');