Я хочу рассчитать среднеквадратичные смещения (MSD) для некоторых частиц в 2D-пространстве. Из того, что я понимаю, MSD является мерой смещений для каждой частицы по траектории: я использую определение, что <(∆r (∆t)) ^ 2> = 1 / N ∑r_i ^ 2 (∆t ) где N - количество частиц.
Смещение рассчитывается как
x_1 = x (t_1), x_2 = x (t_1 + ∆t), ∆x_1 (∆t) = x_2 - x_1
y_1 = y (t_1), y_2 = y (t_1 + ∆t), ∆y_1 (∆t) = y_2 - y_1
...
x_i = x (t_i), x_i + 1 = x (t_i + ∆t), ∆x_i (∆t) = x_i + 1 - x_i
y_i = y (t_i), y_i + 1 = y (t_i + ∆t), ∆y_i (∆t) = y_i + 1 - y_i
Квадратное смещение (∆r) ^ 2 является суммой смещений в каждом измерении. Тогда среднее значение берется.
Как мне это реализовать? Я попробовал следующее, но, как уже отмечали другие, это не правильно.
PROGRAM CALC
IMPLICIT NONE
INTEGER :: J,N,T,NPARTICLES,NSTEPS
REAL(8) :: SUM,DX,DY
REAL(8),ALLOCATABLE :: X(:,:),Y(:,:)
REAL(8),ALLOCATABLE :: MSD(:)
! INPUT
NSTEPS = 101
NPARTICLES = 500
ALLOCATE ( X(NPARTICLES,0:NSTEPS-1) )
ALLOCATE ( Y(NPARTICLES,0:NSTEPS-1) )
ALLOCATE ( MSD(0:NSTEPS-1) )
X = 0.0D0
Y = 0.0D0
DX = 0.0D0
DY = 0.0D0
OPEN(UNIT=50,FILE='TRAJECTORY',STATUS='UNKNOWN',ACTION='READ')
DO T = 0,NSTEPS-1
DO J = 1,NPARTICLES
READ(50,*) X(J,T), Y(J,T)
END DO
SUM = 0.0D0
MSD = 0.0D0
DO WHILE (NSTEPS < T)
DO N = 1,NPARTICLES
DX = X(N,T+1) - X(N,T)
DY = Y(N,T+1) - Y(N,T)
SUM = SUM + (DX**2 + DY**2)
END DO
END DO
MSD(T) = SUM / NPARTICLES
END DO
CLOSE(5)
DEALLOCATE(X)
DEALLOCATE(Y)
OPEN(UNIT=60,FILE='msd.dat',STATUS='UNKNOWN')
DO T = 0,NSTEPS-1
WRITE(60,*) T,MSD(T)
END DO
CLOSE(60)
DEALLOCATE(MSD)
END PROGRAM CALC