Расчет среднеквадратичного смещения с помощью fortran

Я хочу рассчитать среднеквадратичные смещения (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

И я был бы обеспокоен тем, что цикл DO WHILE (NSTEPS > T) не имеет обновленных nsteps или t (если только что-то сумасшедшее не происходит с функциями x или y).

francescalus 31.10.2018 14:16

@idadle, пожалуйста, сделайте полный MWE вашего примера, скомпилируйте его с проверкой границ и запустите. Также посмотрите на цикл WHILE, как уже указывалось francescalus.

albert 31.10.2018 15:58

Не используйте номера единиц 5 и 6 для файлов, поскольку они обычно используются для ввода и вывода через терминал. Лучше использовать функцию newunit (не уверен в названии) или хотя бы количество единиц> 10 (дополнительную информацию ищите на форумах)

albert 31.10.2018 16:20

@albert Спасибо, что указали на это. Сейчас поправлю.

idladle 31.10.2018 16:22

@HighPerformanceMark Спасибо. Я только что отредактировал свой вопрос и код. Что касается смещения, ∆x_i (∆t) = x_i + 1 - x_i, я не уверен, как его рассчитать без петли DO WHILE. Когда я изменяю условие на чтение DO WHILE (NSTEPS < T), оно должно выполняться сейчас, но значения x и y в самом конце не будут иметь значений. Я застрял.

idladle 31.10.2018 16:41

NSTEPS имеет значение 101, а T изменяется от 0 до NSTEPS-1, поэтому от 0 до 100, таким образом, условие NSTEPS <T всегда ложно, и, следовательно, значения DX и DY никогда не вычисляются / сумма не обновляется и т. д.

albert 31.10.2018 16:45

Вы уже запускали проверку границ? (Сейчас это не имеет смысла из-за неправильного цикла while).

albert 31.10.2018 16:49

Отойди от клавиатуры! Отойдите как минимум на 3 метра (если в США ваш Конгресс санкционировал использование метрической системы в 1866 году, поэтому, пожалуйста, не просите меня переводить ее в футы или фарлонги). А теперь остановись и подумай. Подумайте, зачем вы вообще используете цикл do while для перебора структур данных фиксированного размера. Мы знаем, что вы знаете, как использовать «обычный» цикл do, в котором используется индексная переменная. Считать. Не меняйте символы в программе как наугад; Я вижу, как вы меняете > на < и скрещиваете пальцы; что не вышло хорошо не так ли?

High Performance Mark 31.10.2018 16:51

@HighPerformanceMark Спасибо. Было бы неплохо иметь такой оператор, как DO T = 0,NSTEPS-2, вместо цикла DO WHILE? Кажется, это имеет больше смысла в том, что вы говорите о структурах данных фиксированного размера.

idladle 31.10.2018 17:16

Хорошо, теперь, когда вы имеют об этом думаете, вернитесь к клавиатуре, написанию кода и проверьте свои хорошо продуманные идеи. Затем, если что-то не так, задайте вопрос (или отредактируйте его). Но в целом, помните, что я не так заинтересован в вашем успехе, как вы, ваш последний комментарий говорит о том, что вы делаете успехи.

High Performance Mark 31.10.2018 17:46

Помимо логических ошибок (цикл на самом деле не зацикливается), кажется, что ваш код вычисляет среднее значение квадрат смещения за один шаг, усредненное от 0 до T вместо смещения для временного интервала T.

Pierre de Buyl 31.10.2018 22:43
1
11
467
0

Другие вопросы по теме