| 5 | 1/1 | 返回列表 |
| 查看: 4281 | 回復(fù): 8 | |||||
| 當(dāng)前只顯示滿足指定條件的回帖,點(diǎn)擊這里查看本話題的所有回帖 | |||||
[求助]
從LAMMPS得到的原子速度,如何計(jì)算速度自相關(guān)函數(shù)(VAF)?求助 已有3人參與
|
|||||
| 比如說我的系統(tǒng)里有4000個(gè)原子,在某個(gè)狀態(tài)下跑了100步,每一步的時(shí)間間隔是確定的(2fs),每一步時(shí)每個(gè)原子的三維速度都會(huì)輸出。我想計(jì)算這個(gè)系統(tǒng)的速度自相關(guān)函數(shù)VAF,計(jì)算公式是<v(0)·v(t)> 。我不太明白這個(gè)公式的含義是計(jì)算什么,用FORTRAN或MATLAB應(yīng)該如何編程實(shí)現(xiàn)呢?緊急求助,不勝感激。 |
專家顧問 (著名寫手)
![]() |
專家經(jīng)驗(yàn): +218 |
|
% 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)'); |
新蟲 (著名寫手)
銀蟲 (初入文壇)
專家顧問 (著名寫手)
![]() |
專家經(jīng)驗(yàn): +218 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 344求調(diào)劑 +5 | knight344 2026-03-16 | 6/300 |
|
|---|---|---|---|---|
|
[考研] 070300化學(xué)319求調(diào)劑 +6 | 錦鯉0909 2026-03-17 | 6/300 |
|
|
[考研] 302求調(diào)劑 +10 | 呼呼呼。。。。 2026-03-17 | 10/500 |
|
|
[考研] 303求調(diào)劑 +4 | 睿08 2026-03-17 | 6/300 |
|
|
[考研] 278求調(diào)劑 +5 | 煙火先于春 2026-03-17 | 5/250 |
|
|
[考研] 328求調(diào)劑,英語六級(jí)551,有科研經(jīng)歷 +3 | 生物工程調(diào)劑 2026-03-16 | 8/400 |
|
|
[考研] 290求調(diào)劑 +3 | p asserby. 2026-03-15 | 4/200 |
|
|
[考研] 26考研求調(diào)劑 +6 | 丶宏Sir 2026-03-13 | 6/300 |
|
|
[考研] 285化工學(xué)碩求調(diào)劑(081700) +9 | 柴郡貓_ 2026-03-12 | 9/450 |
|
|
[考研] 278求調(diào)劑 +3 | Yy7400 2026-03-13 | 3/150 |
|
|
[基金申請(qǐng)] 國(guó)自科面上基金字體 +6 | iwuli 2026-03-12 | 7/350 |
|
|
[考研] 304求調(diào)劑 +3 | 曼殊2266 2026-03-14 | 3/150 |
|
|
[考研] 070305求調(diào)劑 +3 | mlpqaz03 2026-03-14 | 4/200 |
|
|
[考研] 289求調(diào)劑 +4 | 這么名字咋樣 2026-03-14 | 6/300 |
|
|
[基金申請(qǐng)]
有必要更換申報(bào)口嗎
20+3
|
fannyamoy 2026-03-11 | 3/150 |
|
|
[考研] [0860]321分求調(diào)劑,ab區(qū)皆可 +4 | 寶貴熱 2026-03-13 | 4/200 |
|
|
[考研] 304求調(diào)劑 +7 | 7712b 2026-03-13 | 7/350 |
|
|
[考研] 求b區(qū)學(xué)校調(diào)劑 +3 | 周56 2026-03-11 | 3/150 |
|
|
[考研] 085600材料與化工 309分請(qǐng)求調(diào)劑 +7 | dtdxzxx 2026-03-12 | 8/400 |
|
|
[論文投稿]
投稿問題
5+4
|
星光燦爛xt 2026-03-12 | 6/300 |
|