Кристаллография, 2023, T. 68, № 1, стр. 100-104

ВВЕДЕНИЕ

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

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

МОДИФИЦИРОВАННЫЙ ФУНКЦИОНАЛ ДЛЯ ОПТИМАЛЬНОГО ПРИБЛИЖЕНИЯ ФОРМЫ КРИВОЙ

Одним из важнейших первичных результатов исследования физического явления или процесса является какая-либо экспериментально измеряемая функциональная зависимость, представляемая, как правило, графически в виде кривой. В качестве исходных данных в рефлектометрических исследованиях выступает кривая угловой зависимости интенсивности зеркального рассеяния зондирующего излучения от исследуемого объекта – рефлектометрическая кривая. Обратная задача заключается в нахождении профиля электронной плотности вещества на основе экспериментальной кривой рассеяния [1]. Объектом исследования служит пленка из нескольких слоев жидкокристалличесих G-2 блочных дендримеров со смешанными алифатическими и мезогенными концевыми группами [2]. Рентгеновская рефлектометрическая кривая (зависимость интенсивности зеркального рассеяния от угла скольжения q_z) для G2-дендримера приведена на рис. 1. Зависимость имеет довольно сложную форму, определяемую как внутренней структурой пленки, так и особенностями ее взаимодействия с зондирующим излучением. На кривой рассеяния выделяются начальная угловая область полного внешнего отражения, характерные брэгговские пики, возникающие из-за блочно-периодического строения жидкокристаллической пленки, и осцилляции Киссига, определяемые общей толщиной [3]. При этом наблюдается общее монотонное френелевское падение интенсивности отражения с увеличением угла скольжения, которое составляет приблизительно 6–7 порядков величины [1].

Решение обратной задачи рефлектометрии сводится к поиску такого распределения электронов в поперечном сечении пленки, которое может вызвать распределение интенсивности рассеяния отраженного от пленки излучения, соответствующее экспериментально измеренной кривой. Существует большое количество как модельных, так и безмодельных методов и подходов к решению обратной задачи [1]. В настоящей работе использовался распространенный метод подгонки экспериментальной кривой рассеяния модельной кривой, рассчитанной для ступенчатой модели профиля плотности. Базовым методом расчета интенсивности служил рекуррентный метод Парратта [4].

Процедура нахождения электронного профиля плотности пленки сводится к минимизации определенного функционала (целевой функции), зависящего от параметров конструируемой модели и экспериментально наблюдаемой интенсивности рассеяния.

Наиболее часто используемой целевой функцией в задачах настройки математических моделей на реальные данные является функционал, представляющий собой сумму квадратов отклонений (невязок) рассчитанной величины от экспериментально наблюдаемой. Такой функционал представляет собой критерий качества решения оптимизационной задачи методом наименьших квадратов [5]:

Здесь $I_{i}^{{{\text{exp}}}}$ – экспериментально наблюдаемая величина интенсивности рассеяния, $I_{i}^{{{\text{calc}}}}$ – интенсивность, рассчитанная методом Парратта [4], N – число экспериментальных точек.

Таким образом, находится минимальное значение функции (1) (которое в данном случае равно нулю). Параметры модели пленки, приводящие к значениям {$I_{i}^{{{\text{exp}}}}$}, и будут решением обратной задачи рефлектометрии.

На практике весьма проблематично получить устойчивое адекватное решение, минимизируя функционал (1). Это обусловлено рядом причин. Одна из них – это указанная выше разница в величинах интенсивности начального и конечного угловых диапазонов. Эта разница составляет 6–7 порядков величины, что практически исключает возможность учета вкладов в сумму (1) интенсивностей в области больших углов. Другая важная причина – это сложнейшая форма кривой интенсивности рассеяния (рис. 1), на вид которой дополнительно влияют зашумленность данных и ошибки измерения. Первую сложность (масштаб величин) можно частично учесть демпфированием – формальной заменой величин интенсивностей их логарифмами:

Вопрос подгонки кривой сложной формы остается проблемой при использовании функционалов (1), (2). Для решения этой проблемы предложен функционал специального вида, который реализован в пакете BARD (Basic Analisys of X-Ray Diffraction) [6] для использования в прикладных задачах минимизации при исследовании веществ с помощью методов рентгеновской дифракции. Обозначим через r_i отношение расчетной и измеренной интенсивности в i-й точке и введем следующий функционал F_BARD:

В функционале F_BARD (3) предусмотрены большие штрафы в случае, если вычисляемые значения существенно отличаются от измеряемых величин, т.е. подгоняемая кривая отклоняется от формы оригинальной кривой. Это условие является достаточно жестким, что заставляет даже минимизаторы, основывающиеся на стохастическом методе выбора текущей точки, принимать “правильные” решения. Также отметим, что введение относительной величины r_i (3.1) позволяет избежать описанной выше несоразмерности вкладов величин невязок, вычисленных в различных угловых областях, тем самым приближая кривую целиком к искомой зависимости во всем исследуемом диапазоне.

Рассмотрим применение модифицированного функционала в задаче нахождения электронного профиля плотности пленки по реальным экспериментальным данным. Источником служили рефлектометрические данные пленки блочного G2-дендримера [2]. На рис. 1 приведены результаты подгонки методом отжига рефлектометрической кривой рассеяния от этой пленки с использованием двух функционалов (1) и (3). Видно, что решение на основе модифицированного функционала F_BARD (3) дает результат, наилучшим образом приближенный к экспериментальной зависимости. И, как следствие, определяемый профиль электронной плотности исследуемой пленки оказывается соответствующим всем ее структурным параметрам. Но этот вопрос и его обсуждение остаются за рамками данной работы.

УНИВЕРСАЛЬНАЯ НОРМИРОВКА ПАРАМЕТРОВ ЗАДАЧИ ПОИСКА ЭКСТРЕМАЛЬНЫХ ЗНАЧЕНИЙ

При нахождении электронного профиля пленки обычно используют его представление в виде совокупности слоев различной плотности и толщины (слоистая модель). При этом каждый такой слой характеризуется набором из четырех параметров: толщиной слоя d, поправками к показателям преломления и поглощения (δ, β), связанными с электронной плотностью вещества, и параметром неидеальности (размытия) границы слоя, называемым шероховатостью слоя σ [4]. Параметры d и σ имеют размерность длины, и их значения варьируются от единиц до нескольких тысяч ангстрем для так называемых “толстых” пленок [7]. Безразмерные поправки к показателю преломления n – δ (вещественная часть) и β (мнимая) – связаны с микроскопическими свойствами вещества:

и в рентгеновском диапазоне излучения их величины имеют значения ~10^–6–10^–7 единиц. Таким образом, основные параметры задачи нахождения профиля электронной плотности вещества (обратной задачи рефлектометрии) различаются по абсолютной величине на 6–9 порядков. Из-за такой разницы величин параметров область поиска параметрического N-мерного пространства значений (гиперкуба) представляет собой вытянутые игловидные параллелепипеды, что крайне затрудняет поиск в них экстремумов. Это серьезная проблема большинства современных как локальных, так и глобальных минимизаторов [8]. Например, один из эффективнейших минимизаторов NL2SOL требует задания равномерного “шага” по всем исследуемым параметрам задачи [9]. Чтобы устранить эту проблему, в [10] был разработан специальный метод прямого поиска, не зависящий от масштаба параметров. Также иногда возможен переход к другим физическим единицам измерения, что, вообще говоря, не является универсальным методом.

В данной работе предложен простой и в то же время универсальный способ нормировки переменных при нахождении экстремума для величин различного числового масштаба. В общем виде для решения задачи необходимо найти вектор значений X_P в P-мерном пространстве [X^lo, X^up]^P (lo и up – нижний и верхний пределы), минимизирующий функцию Ф_P(X). Отметим, что простая нормировка параметров на значения своих верхних пределов $X_{P}^{{up}}$ в виде Z_P = X_P/$X_{P}^{{up}}$ (формула обратного пересчета X_P = $X_{P}^{{up}}$Z_P) не приводит к однородной области поиска для всех параметров, и P-мерный параллелепипед остается таковым при этом преобразовании. Это является важным моментом для минимизаторов, в которых в качестве внутренней процедуры поиска используется поиск на сфере [8]. Предложим следующую универсальную формулу пересчета параметров:

При этом область поиска представляет собой единичный P-мерный гиперкуб (значения всех новых переменных Z_P лежат в диапазоне [0–1.0]).

Такого рода нормировка успешно применена для определения по данным дифракции скользящего падения структурной упаковки слоя органических молекул на водной субфазе. Молекулы в форме пластин с размерами сторон в пределах 9 ≤ ≤ a = b ≤ 11 Å и существенно меньшей толщиной c [11] имели тенденцию к “слипанию” по плоскости основания с образованием прямых стопок наподобие кубиков. В слое могли образоваться островки с кристаллической упаковкой, что проявлялось в виде соответствующих пиков на дифракционной картине (рис. 2). В таком случае при поиске структурных параметров слоя учитывались не только различия в числовом масштабе параметров (размеры элементарной ячейки или обратной решетки), но и оценивались такие характеристики дифракционных рефлексов, как индексы Миллера – параметры, принимающие дискретные целочисленные значения из ограниченного диапазона [–10…10]. Поиск структурных параметров для предложенных моделей осуществлялся путем минимизации невязок между координатами дифракционных рефлексов, присутствующих на экспериментальной кривой рассеяния, и аналогичными расчетными значениями, отвечающими целочисленным наборам hkl индексов Миллера, для модельных ячеек при варьировании их параметров в заданных пределах. Таким образом, процедура позволяет определять параметры возможных вариантов структурной организации молекулярных агрегатов в слое даже в случае, когда исследуемое вещество заранее неизвестно [11].

ЗАКЛЮЧЕНИЕ

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

Все расчеты были проведены с помощью аналитического программного комплекса BARD, в рамках которого реализованы предложенные алгоритмы.

Работа выполнена при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию ФНИЦ “Кристаллография и фотоника” РАН.