Фортран 90

PROGRAM dyn_array IMPLICIT NONE INTEGER SIZE
REAL, ALLOCATABLE, DIMENSION(:) ::array WRITE(*, *) 'SIZE?' 
READ(*, *) SIZE
IF(SIZE > 0)  ALLOCATE(array(SIZE))
…
IF(ALLOCATED(array)) DEALLOCATE(array) END PROGRAM dyn_array

PROGRAM NEWTON 
IMPLICIT NONE 
REAL (8) :: X, DX, F, DF 
X = 3.3   !  НАЧАЛЬНОЕ ПРИБЛИЖЕНИЕ DO  !  НЬЮТОНОВСКИЕ ИТЕРАЦИИ
DX = F(X) / DF(X) ! ВЫЧИСЛЕНИЕ ШАГА 
X = X - DX ! ВЫЧИСЛЕНИЕ ОЧЕРЕДНОГО ПРИБЛИЖЕНИЯ 
IF(DX <= SPACING(X)) EXIT ! ЦИКЛ ЗАВЕРШАЕТСЯ,   КОГДА 
!  ШАГ МЕНЬШЕ РАССТОЯНИЯ МЕЖДУ ДВУМЯ ПОСЛЕДОВАТЕЛЬНЫМИ !  ВЕЩЕСТВЕННЫМИ ЗНАЧЕНИЯМИ 
END DO
PRINT *, X !  ВЫВОД ЗНАЧЕНИЯ КОРНЯ
PRINT *, F(X) !  ВЫВОД ЗНАЧЕНИЯ ФУНКЦИИ В ТОЧКЕ X PRINT *, DF(X) !  ВЫВОД ЗНАЧЕНИЯ ПРОИЗВОДНОЙ 
!  ФУНКЦИИ В ТОЧКЕ X
END PROGRAM NEWTON

REAL(8) FUNCTION F(X) IMPLICIT NONE 
REAL(8) :: X
F = SIN(X)
RETURN END

Перечень операторов языка

Ниже дается перечень некоторых конструкций и операторов Фортрана. Необязательные элементы заключены в квадратные скобки []. Если символ пробела не окружен квадратными скобками, то он является обязательным.

Операторы программных компонент

Оператор	Описание
PROGRAM имя_программы	Оператор заголовка главной программы
MODULE имя_модуля	Оператор заголовка модуля
END[ MODULE[имя_модуля]]	Оператор завершения модуля
USE имя_модуля[, ONLY only_список]	Оператор подключения модуля
[RECURSIVE ]SUBROUTINE имя_подпрограммы [([список_формальных_параметров] ) ]	Оператор заголовка подпрограммы-процедуры
[тип][RECURSIVE ]FUNCTION имя_функции ([список_формальных_параметров]) [ RESULT(имя_результата)]	Оператор заголовка подпрограммы-функции
INTERFACE[ родовое_описание]	Оператор заголовка интерфейса
END[ ]INTERFACE	Оператор завершения интерфейса
CONTAINS	Оператор содержания
ENTRY
BLOCK[ ]DATA[ имя_блока_данных]	Оператор заголовка блока данных
END [ ]BLOCK[ ]DATA[имя_блока_данных]	Оператор завершения блока данных

Операторы описания и инициализации данных

Оператор	Описание
MODULE PROCEDURE список_имен_модульных_процедур	Оператор декларации модульных процедур
тип [ [, атрибут] [, атрибут...]....::] список_объектов где тип выбирается из списка: INTEGER[(KIND=]параметр_разновидности_типа)] REAL[(KIND=]параметр_разновидности_типа)] LOGICAL[(KIND=]параметр_разновидности_типа)] COMPLEX[(KIND=]параметр_разновидности_типа)] CHARACTER[список_параметров_типа] DOUBLE[ ]PRECISION] TYPE( имя_типа) и атрибуты - совместимая комбинация из следующих значений: PARAMETER, PUBLIC, PRIVATE, POINTER, TARGET, ALLOCATABLE, DIMENSION(список_экстентов) , INTENT(параметр_входа/выхода), EXTERNAL, INTRINSIC, OPTIONAL, SAVE	Оператор описания
TYPE[, атрибут_доступа ::] имя_производного_типа атрибут_доступа - PUBLIC или PRIVATE	Оператор определения производного типа, заголовок
END[ ]TYPE[ имя_типа]	Оператор определения производного типа: завершение
IMPLICIT список, где список - это тип (список-букв) [, тип (список_букв) ] ... или NONE	Оператор определения правил неявной типизации
ALLOCATABLE [::] имя_массива[( список_экстентов)][, имя_массива [ (список_экстентов) ]...]	Оператор назначения атрибута ALLOCATABLE
DIMENSION имя_массива( список_экстентов)[, имя_массива (список_экстентов)...]	Оператор спецификации массивов
PARAMETER (список_определений_именованных_констант)	Оператор определения именованных констант
EXTERNAL список_внешних_имен	Оператор назначения атрибута EXTERNAL
INTRINSIC список_встроенных_имен	Оператор назначения атрибута INTRINSIC
INTENT(параметр_входа/выхода) список_формальных_параметров	Оператор назначения атрибута INTENT
OPTIONAL список_формальных_параметров	Оператор назначения атрибута OPTIONAL
SAVE[[::] список_сохраняемых_объектов]	Оператор назначения атрибута SAVE
COMMON /[имя_общего_блока]/список_переменных[, /имя_общего_блока/список_переменных.... ]	Оператор описания блока общей памяти
DATA список_объектов/список_значений/[, список_объектов/список_значений/...]	Оператор инициализации объектов
FORMAT([список_дескрипторов])*	Оператор спецификации формата преобразования данных при операциях ввода/вывода

Операторы передачи управления

Оператор	Описание
END[ PROGRAM[ имя_программы]]	Оператор завершения главной программы
END[ подпрограмма[имя_подпрограммы]] где подпрограмма - SUBROUTINE или FUNCTION	Оператор завершения подпрограммы
CALL имя_подпрограммы[(список_фактических_параметров) ]	Оператор передачи управления подпрограмме-процедуре
RETURN	Оператор возврата управления вызывающему компоненту программы
STOP[ сообщение]	Оператор останова программы

Операторы присваивания

Оператор	Описание
переменная = выражение	Оператор присваивания для скалярных и массивоподобных объектов
ссылка => адресат	Оператор прикрепления ссылки к адресату

Управляющие конструкции и условные операторы

Оператор	Описание
IF(скалярное_логическое_выражение) исполняемый_оператор	Условный оператор
WHERE( логическое_выражение_массив) массив_переменная = выражение_массив	Условное присваивание для массивов
[ if_имя : ] IF( скалярное_логическое_выражение) THEN ELSE[[ ]IF(скалярное_логическое_выражение) THEN[ if_имя] END[ ]IF[ if_имя]	Конструкция ветвления IF_THEN_ELSE
WHERE( логическое_выражение_массив) ELSEWHERE END[ ]WHERE	Конструкция ветвления в присваивании массивов
[sе1есt_имя:] SELECT [ ] CASE ( скалярное_выражение) CASE ( список_возможных_значений)[ sе1есt_имя] CASE DEFAULT [ sе1еct_имя] END [ ] SELECT [ sе1еct_имя]	Конструкция SELECT
GO[ ]TO метка	Оператор перехода
[do_имя:] DO[ метка] переменная = скалярное_целое_выражение1, скалярное_целое_выражение2[, скалярное_целое_выражение3]	Заголовок оператора цикла
[do_имя:] DO[ метка] [,] WHILE(скалярное_лог_выражение)	Заголовок оператора цикла в альтернативной форме
CYCLE [ do_имя]	Оператор перехода к оператору завершения цикла с именем
EXIT[ do_имя]	Оператор выхода из цикла с именем do_имя
CONTINUE	Оператор перехода к следующему шагу цикла
END[ ]DO[ do_имя]	Оператор завершения цикла с именем do_имя

Операторы управления динамической памятью

Оператор	Описание
ALLOCATE( список_выделяемых_объектов[, STAT=статус])	Оператор выделения динамической памяти для выделяемых объектов
DEALLOCATE( список_выделенных_объектов[, STAT=статус])	Оператор освобождения динамической памяти от выделенных объектов

Операторы ввода-вывода

Оператор	Описание
READ( список_управления_вводом) [ список_ввода] READ формат[, список_ввода]	Оператор чтения данных
WRITE(список_управления_выводом) [ список_вывода]	Оператор записи данных
PRINT формат[, список_вывода]	Оператор вывода данных на устройство стандартного вывода
OPEN( список_спецификаций)	Оператор соединения файла с логическим устройством ввода-вывода
CLOSE(список_спецификаций)	Оператор закрытия файла

Встроенные подпрограммы

Имя подпрограммы	Возвращаемое значение или производимое действие
ABS(A)	Абсолютная величина
ACHAR(I)	I-ый символ сортирующей последовательности ASCII
ACOS(X)	Арккосинус в радианах
AIMAG(Z)	Мнимая часть комплексного числа
AINT(A[, KIND])	Усечение до целого
ALLOCATED(ARRAY)	Проверка выделенности динамического массива
ANINT(A [, KIND])	Ближайшее целое
ASIN(X)	Арксинус в радианах
ATAN(X)	Арктангенс в радианах
ATAN2(Y, X)	Аргумент комплексного числа
CALL DATE_AND_TIME([DATE][,TIME] [,ZONE] [,VALUES])	Считывание времени и даты с часов реального времени
CALL RANDOM_NUMBER(HARVEST)	Случайное число в интервале [0,1)
CALL RANDOM_SEED([SIZE] [,PUT][,GET])	Получение/назначение затравочного массива
CALL SYSTEM_CLOCK([COUNT] [,COUNT_RATE] [,COUNT_MAX])	Целочисленный отсчет часов реального времени
CEILING(A)	Наименьшее целое не меньшее аргумента
CHAR(I [, KIND]) и ICHAR(C)	I-ый символ в сортирующей последовательности процессора
CMPLX(X [, Y] [, KIND])	Конструктор комплексного числа
CONJG(Z)	Комплексное сопряжение
COS(X)	Косинус
COSH(X)	Гиперболический косинус
COUNT(MASK [, DIM])	Число элементов массива, имеющих значение "истина"
CSHIFT(ARRAY, SHIFT [, DIM])	Циклический сдвиг элементов массива
DIGITS(X)	Число значащих цифр в модели представления аргумента
DOT_PRODUCT(VECTOR_A, VECTOR_B)	Скалярное произведение векторов
DPROD(X, Y)	Умножение с двойной точностью
EOSHIFT(ARRAY, SHIFT [, BOUNDARY] [, DIM] )	Вытесняющий сдвиг элементов массива
EPSILON(X)	Наименьшее число в модели представления аргумента, не пренебрежимое по сравнению с единицей
EXP(X)	Экспонента
EXPONENT(X)	Степенная часть в модели представления аргумента
FLOOR(A)	Наименьшее целое, не превышающее аргумент
FRACTION(X)	Дробная часть в модели представления аргумента
HUGE(X)	Наибольшее число в модели представления аргумента
IACHAR(С)	Индекс символа-аргумента в сортирующей последовательности ASCII
IAND(I, J)	Побитовое логическое И
IBCLR(I, POS)	Установка нуля в заданном бите
IBITS(I, POS, LEN)	Извлечение цепочки битов
IBSET(I, POS)	Установка единицы в заданном бите
ICHAR(C)	Индекс символа-аргумента в сортирующей последовательности процессора
IEOR(I, J)	Побитовое ИСКЛЮЧАЮЩЕЕ ИЛИ
INDEX(STRING, SUBSTRING [,BACK])	Индекс начала подстроки в строке
INT(A [, KIND])	Приведение к целому типу
IOR(I, J)	Побитовое логическое ИЛИ
ISHIFT(I, SHIFT)	Логический сдвиг битов
ISHIFTC(I, SHIFT [, SIZE])	Логический циклический сдвиг части битов вправо
KIND(X)	Разновидность типа аргумента
LBOUND(ARRAY [, DIM])	Нижняя граница массива
LEN(S)	Длина текстового аргумента
LEN_TRIM(STRING)	Длина строки без учета конечных пробелов
LOG(X)	Натуральный логарифм
LOG10(X)	Десятичный логарифм
LOGICAL(L [, KIND])	Приведение к заданной разновидности логического типа
MATMUL(MATRIX_A, MATRIX_B)	Перемножение матриц
MAX (A1, A2 [, A3, ])	Выбор максимального аргумента
MAXLOC(ARRAY [, MASK])	Индекс наибольшего элемента массива
MAXVAL(ARRAY [, DIM ] [,MASK])	Наибольший элемент массива
MIN (A1, A2 [, A3, ])	Выбор минимального аргумента
MINLOC(ARRAY [, MASK])	Индекс наименьшего элемента массива
MINVAL(ARRAY [, DIM ] [,MASK])	Наименьший элемент массива
MOD(A, P)	Остаток от деления по модулю
MODULO(A, P)	Деление по модулю
NINT(A [, KIND])	Ближайшее целое
NOT(I)	Побитовое логическое дополнение
PRECISION(X)	Десятичная точность в модели представления аргумента
PRESENT(A)	Проверка присутствия необязательного аргумента
PRODUCT(ARRAY [, DIM] [,MASK])	Произведение элементов массива
REAL(A [, KIND])	Приведение к вещественному типу
RESHAPE(SOURCE, SHAPE [,PAD] [, ORDER])	Изменение формы массива
SCAN(STRING,SET[,BACK])	Индекс крайнего символа строки STRING в наборе SET
SELECTED_INT_KIND(R)	Параметр разновидности целого типа для заданного степенного диапазона
SELECTED_REAL_KIND([P][,R])	Параметр разновидности вещественного типа для заданной точности и/или заданного степенного диапазона
SHAPE(SOURCE)	Форма аргумента
SIGN(A, B)	Абсолютное значение A со знаком B
SIN(X)	Синус
SINH(X)	Гиперболический синус
SIZE(ARRAY[,DIM])	Размер массива
SQRT(X)	Квадратный корень
SUM(ARRAY[,DIM][,MASK])	Сумма элементов массива
TAN(X)	Тангенс
TANH(X)	Гиперболический тангенс
TINY(X)	Наименьшее положительное число в модели представления аргумента
TRANSPOSE(MATRIX)	Транспонирование матрицы
UBOUND(ARRAY [, DIM])	Верхняя граница массива

Лабораторная работа 3.1 Оптимизация и распараллеливание вычислительной программы на примере моделирования системы взаимодействующих частиц методом молекулярной динамики

Метод молекулярной динамики

Метод молекулярной динамики является одним из основных методов моделирования динамики систем, состоящих из взаимодействующих между собой частиц. Он применяется в молекулярной физике, астрофизике и других науках. Интерпретация понятия "частица" зависит от предметной области. В молекулярной физике это могут быть атом или молекула, а в астрофизике "частицами" являются звезды или планеты, взаимодействующие с другими астрономическими объектами посредством гравитационного взаимодействия. Динамика систем таких частиц часто описывается уравнениями классической механики.

Будем рассматривать систему N частиц, динамика которых описывается уравнением второго закона Ньютона:

$$F_i=m_i a_i$$

Здесь $$F_i$$ - равнодействующая всех сил, действующих на $$i$$ -ю частицу, $$m_i$$ - ее масса и $$а_i$$ - ускорение, с которым частица двигается под действием силы. Это обыкновенное дифференциальное уравнение второго порядка, которое можно записать в виде:

$$F_i=m_i\frac{d^2 r_i}{dt^2}$$

где $$r_i$$ радиус-вектор $$i$$ -й частицы.

Равнодействующая сил, является суммой парных взаимодействий $$i$$ -й частицы со всеми остальными частицами:

$$F_i=F(r_i)=\sum\limits_{i\ne j=1}^N F(r_i,r_j)$$

или

$$F_i=\sum\limits_{j=1}^{i-1}F(r_i,r_j) + \sum\limits_{j=i+1}^N F(r_i,r_j)$$

Сила взаимодействия связана с потенциалом взаимодействия. В расчетах молекулярных систем часто используют потенциал Леннарда-Джонса:

$$V(r)=4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]$$

Здесь $$r =\mid r\mid $$, а $$\epsilon$$ и $$\sigma$$ - параметры потенциала. Первый определяет интенсивность взаимодействия, а второй - положение минимума (см. рис. 6.5).

(рис 6.5) Потенциал Леннарда-Джонса

Потенциал Леннарда-Джонса хорошо описывает взаимодействие между атомами Аргона. Для аргона $$m\approx 6.69\times 10^{-23}\hat e \tilde a,\:\:\: \sigma\approx 3.405 \times 10^{-10} \grave i ,\:\:\: \epsilon \approx 1.654\times 10^{-21} \ddot A ae$$. Будем считать, что массы всех частиц одинаковы, система единиц выбрана таким образом, что $$m = 1$$, $$\epsilon=1$$ и $$\sigma=1$$. Тогда сила парного взаимодействия:

$$F(r_i,r_j)=\frac{24(r_i-r_j)}{{\mid r_i-r_j\mid }^2} \left[\frac{2}{{\mid r_i-r_j\mid }^{12}} - \frac{1}{{\mid r_i-r_j\mid }^6} \right]$$

Интерес могут представлять траектории частиц или термодинамические характеристики системы, например, зависимость температуры от времени:

$$T(t)=\frac{1}{3Nk_B}\sum\limits_{i=1}^N {\left \mid \frac{dr_i(t)}{dt} \right \mid }^2$$

где $$k_B\approx 1.38\times 10^{-23}\frac{\ddot A ae}{\hat E}$$ -постоянная Больцмана.

Для вычисления траекторий в методе молекулярной динамики используются различные алгоритмы. Они основаны на замене непрерывного времени дискретным набором значений $$t\to t^k=t^0+k\Delta t, \:\:\: k=1,2,\ldots$$.. Одним из алгоритмов является метод Верле:

$$r_i^{k+1}=2r_i^k-r_i^{k-1}+F(r_i^k)\frac{\Delta t^2}{m_i}$$

Здесь $$r_i^k=r_i(t^k)$$. Наиболее трудоемким является вычисление силы, действующей на частицу.

Для оптимизации вычисления сил иногда вводится радиус обрезания - некоторый фиксированный (достаточно большой) радиус, такой, что пересчет сил, действующих между частицами, оказавшимися внутри этого радиуса, производится на каждом шаге. Пересчет сил, действующих со стороны более удаленных частиц, производится через несколько шагов. Такая стратегия позволяет уменьшить количество трудоемких операций расчета парных сил. Номера "ближайших" частиц могут храниться в специальном массиве. прТакой массив называется списком.Обычно используются два массива целого типа. Каждый элемент первого массива ( list_nb ) содержит количество ближайших соседей соответствующей частицы:

k0 = 0 
j0 = 0 
do i = 1,  n 
do j = 1, n
if(abs(x(i) - x(j)) < rcutoff)   then 
j0 = j0 + 1
list_nb(i)  = j end if 
end do
i_nb(i) = j0 - k0 end do

Во втором массиве содержатся номера соседей каждой частицы. Положения хранятся в элементах, начиная с start(i) и заканчивая finish(i):

start(1) = 1
finish(1) = i_nb(1) do i = 2, n
start(i) = finish(i - 1) + 1
finish(i)  = finish(i - 1) + i_nb(i) end do

Лабораторная работа

В заданиях лабораторной работы 3.1 предлагается выполнить оптимизацию вычислительной программы. Первый этап оптимизации - сокращение числа операций в последовательном коде на основе использования физических законов. Затем выполняется распараллеливание программы. Выбор инструмента распараллеливания должен быть аргументирован.

Данная работа не предполагает выполнения всех заданий. Учащемуся может быть предложен набор из нескольких заданий, обязательно включающий номера 1, 2 и 3, а также одно или несколько заданий, связанных с распараллеливанием.

По результатам составляется отчет, содержащий решения с обоснованием выбранного метода, результаты измерения быстродействия программы, а также физические результаты и выводы.

Выполнение данной работы рассчитано на несколько занятий (2-6 академических часов). Имеется полный исходный текст программы на языке Fortran 90.

Задания для практической работы

Задание 1

Получить у преподавателя файл с исходным текстом программы и ознакомиться с реализацией метода молекулярной динамики.

Задание 2

Откомпилировать программу, выполнить моделирование для системы с параметрами, заданными преподавателем. Определить процессорное время, потраченное на выполнение моделирования.

Задание 3

Выполнить оптимизацию последовательного кода программы на основе "физических соображений". Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в задании 2.

Задание 4

Проанализировать последовательный код и выявить участки потенциального параллелизма. Выполнить распараллеливание с помощью OpenMP (если это возможно). Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в задании 3.

Задание 5

Проанализировать последовательный код, выявить участки потенциального параллелизма и предложить способ распараллеливания в модели передачи сообщений. Выполнить распараллеливание с помощью MPI, используя блокирующие двухточечные операции обмена (если это возможно). Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатами, полученными в заданиях 3 и 4.

Задание 6

Проанализировать возможность использования других видов обмена (двухточечные неблокирующие, коллективные). Если предполагается выигрыш в производительности, модифицировать вариант программы с использованием MPI. Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатами, полученными в заданиях 3 , 4 и 5.

Задание 7

Выполнить оптимизацию последовательного кода, полученного при выполнении задания 3, используя введение радиуса обрезания и метод генерации списков частиц. Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в задании 3.

Задание 8

Провести анализ зависимостей при генерации списков частиц. Выполнить распараллеливание, выбрав наиболее подходящий инструмент. Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в задании 7.

Задание 9

Выполнить распараллеливание вычисления сил, действующих на частицы, используя метод списков. Определить процессорное время, потраченное на выполнение моделирования. Сравнить с результатом, полученным в заданиях 4 или 5.

Задание 10

Исследовать степень сбалансированности загрузки процессоров (ядер при использовании многоядерной архитектуры). Если сбалансированность неудовлетворительная, предложить способы ее улучшения.

Задание 11

Исследовать масштабируемость программы. Если масштабируемость неудовлетворительная, предложить способы ее улучшения.

Задание 12

На основании результатов, полученных при выполнении заданий данной лабораторной работы, написать отчет, в котором содержатся выводы об эффективности различных способов оптимизации исходного последовательного кода и трудоемкости реализации этих способов на практике.

Пример 1

Далее приводится листинг программы молекулярной динамики.

program mol_dyn

implicit real(8) (a-h, o-z) 
integer, parameter :: ndim = 1000 
!
! Массивы координат,  проекций скоростей и ускорений частиц 
!
real(8), dimension(1:ndim) :: x, y, vx, vy, ax, ay
!
!  Задание начальной конфигурации частиц
!
call start(x, y, vx, vy, N, Sx, Sy, dt, dt2, nsnap, ntime) 
!
!  Вычисление ускорений частиц
!
call accel(x, y, ax, ay, N, Sx, Sy, zpe)

zpe =0.0d0
!
!  Энергия выводится через nsnap шагов
!
do isnap = 1, nsnap
!
! Внутренний цикл - по шагам по времени
!
do itime = 1, ntime
!
!  Выполняется сдвиг частиц
!
call move(x, y, vx, vy, ax, ay,N, Sx, Sy, dt, dt2, zke, zpe) 
end do
!
! Вывод значений энергии:  кинетической,  потенциальной и полной
!
call output(zke,   zpe,   Sx,   Sy,  dt,  N,  ntime)

end do

end program !
!====================================================
!
subroutine start(x,  y,  vx,  VY,  N,   Sx,   Sy,  dt,  dt2,  nsnap,  ntime)

implicit real(8) (a-h, o-z) integer,  parameter :: ndim = 1000
real(8), dimension(1:ndim) ::  x,  y, vx, vy, ax, ay
!
!  Число частиц 
!
n = 200
!
!  Размер области моделирования
!
sx = 100.0d0
sy = 100.0d0 !
! Шаг по времени
!
dt = 1.0d-11 dt2 = dt**2 !
! Максимальное значение скорости
!
vmax = 10.0d0 !
!   Число строк в таблице вывода и количество шагов по времени между ними
!
nsnap = 200
ntime = 500
!
!  Формируется начальная хаотическая конфигурация
!
do i = 1, N
call rndm(ranx) 
x(i) = sx * ranx 
call rndm(rany) 
y(i) = sy  * rany 
call rndm(ranx) 
call rndm(rany) 
vx(i) = vmax * (2  * ranx - 1) vy(i)  = vmax * (2 * rany - 1) end do

do i = 1, N
vxcum = vxcum + vx(i) vycum = vycum + vy(i) 
end do

vxcum = vxcum / N vycum = vycum / N

do i = 1, N
vx(i) = vx(i) - vxcum vy(i) = vy(i) - vycum 
end do

end subroutine

!
!====================================================
!
subroutine move(x, y, vx, vy, ax, ay, N, Sx, Sy, dt, dt2, zke, zpe)

implicit real(8) (a-h, o-z) 
integer, parameter :: ndim = 1000
real(8), dimension(1:ndim) ::  x,   y,   vx,   vy,   ax,   ay

do i = 1, N
xnew = x(i) + vx(i) * dt + 0.5d0 * ax(i) * dt2 
ynew = y(i) + vy(i) * dt + 0.5d0 * ay(i) * dt2
!
! Учет периодических граничных условий
!
call cellp(xnew,  ynew, vx(i), vy(i), sx, sy) 
x(i) = xnew 
y(i) = ynew
vx(i) = vx(i) + 0.5d0 * ax(i) * dt vy(i) = vy(i) + 0.5d0 * ay(i) * dt 
end do

call accel(x, y, ax,  ay,  N,   Sx,   Sy,   zpe)

do i = 1,  n
vx(i) = vx(i)  + 0.5d0  * dt * ax(i)
vy(i) = vy(i)  + 0.5d0  * dt * ay(i)
zke = zke + vx(i)**2 + vy(i)**2 end do

end subroutine
!
!====================================================
!
subroutine accel(x, Y, ax, ay, N, sx, sy, zpe)

implicit real(8) (a-h, o-z)
integer,  parameter  ::  ndim = 1000
real(8),  dimension(1:ndim) :: x, y, ax, ay

do i = 1,  n
 ax(i) = 0.0d0
ay(i) = 0.0d0 end do

do i = 1,  n 
do j = 1,  n 
 if(i /= j) 
  then 
  dx = x(i) - x(j)
dY = y(i) - y(j)

if(dabs(dx)>0.5 * sx) dx = dx - sign(sx, dx) if(dabs(dy)> 0.5 * sy)  dy = dy - sign(sy,  dy)

!  Вычисление силы,  действующей на j-ю частицу
!
r = dsqrt(dx**2 + dy**2) 
ri = 1.0d0 / r 
ri6 = ri**6
g = 24.0d0 * ri6 * ri * (2.0d0 * ri6 - 1)
force = 0.5d0 * g * ri
pot = 4.0d0  * ri6  * (ri6**2 - 1.0d0)

ax(i) = ax(i) +	force	* dx
ay(i) = ay(i) +	force	* dy
ax(j) = ax(j) -	force	* dx
ay(j) = ay(j) -	force	* dy
zpe = zpe + pot
end if

end do 
end do

end subroutine !
!====================================================
!
subroutine cellp(xnew,  ynew,  vx,  vy,   Sx,   Sy)

implicit real(8)   (a-h,   o-z)
if(xnew < 0) xnew = xnew + sx
if(xnew > sx) xnew = xnew - sx
if(ynew < 0) ynew = ynew + sy
if(ynew > sy)ynew = ynew - sy

end subroutine
!
!====================================================
!
subroutine output(zke, zpe, sx, sy,dt, n, ntime) 
implicit real(8)   (a-h,   o-z)
data iff/0/

if(iff == 0) then iff = 1
write(6,*) '	ke pe tot'
end if

zke = 0.5d0 * zke / ntime zpe = zpe  / ntime 
tot = zke + zpe

write(6, "(6(1x,  e13.6))") zke, zpe,tot

zke = 0.0d0
zpe = 0.0d0

end subroutine 
!
!  Генератор псевдослучайных чисел 
!
subroutine rndm(ran) implicit real*8(a-h,o-z)
data a31/2147483647.d0/,  mul/16807.d0/,  x/268435451.d0/ 
t = dmod(x * mul, a31) 
ra = t
x = t
ran = ra / a31 end subroutine

Приложение. Рекомендации по выполнению работы

Задание 1

Исходный текст последовательного варианта программы снабжен комментариями.

Задание 3

Применение третьего закона Ньютона:

$$F(r_i,r_j)=- F(r_i,r_j)$$

позволяет сократить объем вычислений в цикле примерно в 2 раза.

Следующий прием основан на введении "радиуса обрезания" $$R_{cut}$$ и разбиении всех частиц на две группы по отношению к данной:

частицы, для которых $$\mid r_i-r_j\mid \le R_{cut}$$ ;

частицы, для которых $$\mid r_i-r_j\mid > R_{cut}$$

$$R_{cut}$$ выбирается таким образом, чтобы частицы из первой группы давали основной вклад во взаимодействие. Выбор обычно выполняется эмпирически, то есть, на основе пробных расчетов.

При вычислении силы сумма делится на две подсуммы по группам частиц: $$F_i=\sum\limts_{\mid r_i-r_j \mid \le R_{cut}} F(r_i,r_j) + \sum\limts_{\mid r_i-r_j \mid > R_{cut}} F(r_i,r_j)$$

Если вклад от удаленных частиц пересчитывать реже, чем от частиц, для которых $$\mid r_i-r_j\mid \le R_{cut}$$, трудоемкость вычисления силы уменьшится.