Суперкомпьютерные методы моделирования связанной квантовой системы: кубит-измерительный осциллятор*

Abstract

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

-

нений, мы используем квантовый метод Монте-Карло, который позволяет решать N  2 уравнений для стохастической волновой функции [8]. Это означает, что необходимо многократно решать уравнение Шредингера со случайными параметрами, а наблюдаемые величины определять путем усреднения билинейной комбинации монте-карловской волновой функции по большому числу реализаций (число определяется необходимой точностью). Решение независимых дифференциальных уравнений хорошо параллелизуется, поэтому можно ожидать значительного ускорения процесса вычисления усредненных величин, аналогичным тем, что получают при решении матрицы плотности. В работе описаны принципы распараллеливания и численного моделирования стохастического уравнения Шредингера, а также программный комплекс, позволяющий осуществлять высокопроизводительные параллельные расчеты на GPU. 2. Модель системы "кубит+осциллятор"

В качестве примера измерительного устройства выбран джозефсоновский осциллятор, взаимодействующий с зарядовым кубитом специального вида (так называемый “квантрониум кубит” - quantronium qubit) [5, 6, 9]. Такой кубит работает на основе двух хорошо известных явлений, которые подробно описаны в литературе [10]. Во-первых, это эффект Джозефсона, имеющий место при протекании сверхпроводящего тока через тонкую диэлектрическую прослойку [10]. Во-вторых, используется эффект кулоновской блокады, возникающий вследствие существования конечной энергии, которую нужно затратить, чтобы перенести заряд между контактом и маленьким сверхпроводящим островком [5, 6].

Предполагается, что два джозефсоновских перехода соединены со сверхпроводящим островком, функционирующим в режиме кулоновской блокады, имеют полную емкость CJ , и каждый из них обладает джозефсоновской энергией EJ / 2 [5, 6]. Для управления зарядом на островок подается постоянный потенциал Vg (см. рис. 1). Рис. 1. Схема кубита (quantronium qubit) (центр рисунка), связанного с устройством записи (левая часть рисунка) и взаимодействующего со слабодиссипативным измерительным джозефсоновским осциллятором (правая часть рисунка). Перечеркнутыми квадратиками обозначены джозефсоновские переходы. Островок (кружок, разделяющий перечеркнутые квадраты) является в то же время вторым контактом переходов (см. [5,6]).

Полный гамильтониан системы включает гамильтониан расщепленного куперовского островка

HQQ  4EC  Nˆ  1  CgV (t) 2 

  EJ cos cos ,  2 2e  2 а гамильтониан измерительного джозефсоновского осциллятора имеет вид:

Qˆ 2 2 2C HJ   EJR cos  2e

I (t) , где Nˆ и Qˆ / 2e обобщенные импульсы, сопряженные к фазам  и  , соответственно купе(1) ( 2 ) ровского и джозефсоновского переходов. Управление состояниями кубита осуществляется переменным напряжением V (t)  Vg Vrf (t) cost , а переменный ток I (t) служит для управления измерительным джозефсоновским осциллятором.

Вблизи точки кулоновского вырождения можно использовать базисные состояния | N  и | N  1  ( Nˆ | N  N | N  ) для состояний островка. В этом подпространстве: | N  (1,0) , | N  1  (0,1) , а гамильтониан аппроксимируется выражением

CgVrf (t)  Q2 2 H  2EC  x  EJ cos 2  z  2C  EJR cos  I (t) , ( 3 ) e 2e где  x и z - спиновые матрицы Паули.

При относительно малом токе I (t) джозефсоновский осциллятор работает в слабо нелинейном режиме. Тогда потенциал в выражение ( 3 ) можно разложить по степеням  :

H  2EC CgVerf (t) x  EJ z  2QC2  EJ 1   z 22  EJ 1  4  z 44!  2e2 I (t) , ( 4 ) где   EJ / 4EJR . Вводя операторы рождения a и уничтожения a для параметризации координаты 1/4 1/ 4 и импульса осциллятора

Nˆ  i(a  a ) E R / 2EJ 

C ( 6 ) (7) E R

c , 12 ( Qˆ  2eNˆ , Nˆ  i   ) с учетом  , Nˆ   i , можно гамильтониан ( 4 ) связанной системы "кубит+нелинейный осциллятор" представить в виде составляющих:

H  Hq  Hosc  Hint . ( 5 ) В выражении ( 5 ) первое слагаемое - гамильтониан кубита, который запишем в стандартном виде:

Hq      (t) x  .

q z 2 Можно видеть, что расстояние между уровнями кубита ( q  2EJ ) определяется частотой  q / 2 (~10 ГГц ); а управляющая функция  (t)  4EC e реходы между состояниями кубита  ( z    ,   1). Предполагается, что управление кубитом будет осуществляется Раби-импульсами переменного поля  (t)  Acost определенной длительности, позволяющими приготавливать наперед заданную суперпозицию состояний кубита.

Таким образом, измерение состояний кубита будет осуществляться нелинейным осциллятором, будет индуцировать пеHosc   aa   (a  a)4  f (t)(a  a ) , 0

CgVrf (t) где собственная частота осциллятора 0 

 2EcR 1/ 4 f (t)  f0 cos t , где амплитуда f0   и частота  внешнего возмущения.

2e  EJ  Наконец, последнее слагаемое в ( 5 ) отвечает за взаимодействие кубита с измерительным осциллятором: 2eIc , параметр нелинейности  

C Hint    0aa  1 4  (a  a )4  z , где  - константа связи кубита и осциллятора. Отметим, если управляющее поле отсутствует,  (t)  0 , то оператор взаимодействия ( 8 ) коммутирует с гамильтонианом кубита ( 6 ), поэтому осциллятор производит так называемое «неразрушающее измерение» [5].

Для системы "кубит+осциллятор" уравнение для оператора плотности  в форме Линдблада, в борн-марковском приближении [7] записывается в виде:  t  1 i ( 8 ) , ( 9 )   2

(2a a  aa   aa) где    ( x  i y ) / 2 , а   и  e - соответственно фазовая и энергетическая скорости релаксации кубита,  - параметр релаксации осциллятора. Отметим, что данные параметры выражаются стандартным образом через корреляционными функциями бозонного термостата ( Fz (t) Fz , R(t) R и т.д.) (см. вывод в [7]). 3. Метод квантовых траекторий (квантовый метод Монте-Карло) для системы "кубит+осциллятор"

Уравнение для матрицы плотности ( 9 ) дает усредненную по ансамблю динамику системы. Чтобы проследить динамику кубита и осциллятора в каждом акте измерения (аналог "однократных" измерений) обратимся к эквивалентону языку описания динамики - методу квантовых траекторий или, как его ещё называют, квантовому методу Монте-Карло [8]. Согласно данному подходу эволюция системы, может быть переформулирована как динамика состояний «квантовых траекторий». 3.1. Алгоритм расчетов квантовым методом Монте-Карло для системы "кубит+осциллятор"

Применяя конечно-разностную аппроксимацию для производной по времени в уравнении ( 9 ), перепишем его в виде:  (t  t)  U (t)U   t   z z   e     a a  , (10) где первое слагаемое в ( 2 ) определяется оператором эволюции U  eiHdist / диссипативную динамику системы с гамильтонианом Вигнера-Вайскопфа [7]   

Hdis  Hq  i 4  z z  i 2e     i 2 aa . ( 11 ) Второе же слагаемое в уравнении (10) можно связать со спонтанными переходами (скачками) в системе за счет взаимодействия с бозонным термостатом. Рассмотрим некоторую реализацию и отвечает за чистого состояния  (t)   (t)  (t) . Тогда из (10) можно увидеть, что изменение матрицы плотности  (t  t)   (t  t)  (t  t) обусловлено двумя возможными вкладами: 1) либо с вероятностью 1  P согласно

, 1  P происходит изменение волновой функции (“траектории”) за счет диссипативной динамики (здесь введены следующие обозначения P  t ),       z z   e        aa  ); 2) либо с вероятностью P изменение “траектории” сопровождается одним из скачков:  s (t)   z  (t)   z z  ,  es (t)     (t)      и  s (t)  a  (t)  aa  , ( 13 ) взятыми соответственно с весами 

 p     z z  , pe  e        и p   aa  .

Для нахождения усредненной динамики (аналог того, что получаем в результате решения уравнения ( 11 )) следует получить набор M единичных реализаций, а далее усреднить величину по ним 1 M     i (t)  i (t) . ( 13 )

M i1 После того, как матрица плотности найдена, можно вычислить ожидаемое значение любой наблюдаемой квантовой системы. Например, среднее число квантов осциллятора (средняя энергия): n  i1 n0

M N   n Cni (t)

2 M

,

N 1  i (t)  Cni (t) n .

n0 (12) (14) где Cni (t) - коэффициенты разложения "монте-карловской" функции по базису, образуемому собственными функциями оператора чисел заполнения nˆ n  n n , nˆ  aa : 3.2. Схема распараллеливания квантового метода Монте-Карло

Расчет физических величин может быть организован с использованием GPU ускорителей с применением технологии CUDA [11-12]. Поскольку реализации статистически независимы, отсутствует необходимость в обмене данными между блоками и потоками внутри отдельных блоков, создаваемых на GPU, то есть на каждом процессоре выполняется одна и та же подпрограмма, что позволяет обеспечить параллелизм на уровне данных (SIMT). Однако так как для достижения точности метода МК в 1% необходимо произвести расчет эволюции системы с учетом усреднения по 10000 квантовым траекториям. Поэтому для большей производительности программы эффективно использовать в расчетах кластерную систему, узлы которой содержат GPU. В частности, вычислительные эксперименты проводились на кластере "Лобачевский" Нижегородского государственного университета им. Н.И.Лобачевского [13]. Мы использовали до 64 узлов следующей конфигурации: 2x Intel Xeon E5-2660 CPU (8 cores, 2.2 GHz), 64 GB RAM, NVIDIA Kepler K20Х (2 или 3), OS Linux (CentOS 6.4), Intel C/C++ Compiler, Intel MPI, CUDA 6.4 Toolkit. На каждом узле использовались 2 графических сопроцессора. Распределение вычислительной нагрузки между узлами кластера выполнено с использованием технологии MPI. На каждом вычислительном узле кластера для расчета траекторий используются возможности графических ускорителей. Поскольку расчет отдельных квантовых траекторий является независимым, каждая из них рассчитывается отдельным потоком GPU, при этом один поток GPU может последовательно рассчитать несколько траекторий. Незначительная доля последовательного кода не оказывает существенного влияния на время выполнения.

Для использования всех имеющихся графических ускорителей работа с данными организована следующим образом:

 на управляющем хосте (нулевой процесс) запускается основной модуль программы, который инициализирует данные о параметрах системы, вычисляет начальное состояние, выполняет инициализацию MPI, а также создает структуры необходимых библиотек;  считываются данные о количестве реализаций (квантовых траекториях) и создаются несколько копий структуры с информацией о гамильтониане, которые будут использоваться на отдельных параллельных потоках, исполняемых на каждом из используемых GPU, причем данные разбиваются на пачки (для каждого GPU) таким образом, чтобы на один блок приходилось по 256 потоков, а количество блоков CUDA определяет самостоятельно запрашивая необходимые данные о встроенных GPU-ускорителях;

 с помощью технологии MPI, с исполняемого хоста, создаются параллельные потоки на каждом вычислительном узле, в соответствии с числом GPU устройств, каждый из которых имеет свой уникальный индетификационный номер. Количество исполняемых заданий на каждой видеокарте соответствует числу квантовых траекторий, по которым производится усреднение. Для этого общее число траекторий делится на части: m1, m2, ..., mM, каждая из которых исполняется на своем GPU.

 инициализируется ядро вычислений на каждом из GPU и запускается алгоритм программы, промежуточные результаты накапливаются в памяти каждого из используемых графических адаптеров;

 после завершения расчетов всех траекторий на отдельной GPU происходит усреднение данных (т.е. соответственно по m1, m2, ..., mM);

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

Одним из важных условий эффективности работы программы на GPU является оптимальное использование имеющихся типов памяти: глобальной, разделяемой и памяти констант. В частности для решения данной задачи необходимы следующие наборы данных: информация о гамильтониане (хранятся в регистрах разделяемой памяти), константы (хранятся в кешируемой памяти констант), текущие населённости уровней (хранятся в регистрах разделяемой памяти) и результирующие данные (имеют большой объем, хранятся в глобальной памяти). Регистры являются быстрым типом памяти, но имеют небольшой объем, поэтому использование их для работы с результирующими данными не представляется возможным. Глобальная память, напротив, позволяет хранить большие объемы данных, но является относительно медленной. Однако в данной задаче обращения к результирующим данным массивов, находящимся в глобальной памяти, немного, поэтому хранение в глобальной памяти не оказывает существенного влияния на производительность программы. А также отметим, что если при усовершенствовании архитектуры GPU объем разделяемой памяти будет увеличен, то существует возможность хранить результирующие данные в данной памяти GPU, что позволит ещё более эффективно использовать обращение к памяти и ускорить расчет.

Заметим, что точность квантового метода Монте-Карло зависит от числа реализаций ~ 1 M , где M – число реализаций. В силу наличия стохастических процессов двойная точность для расчётов не требуется, а использование операций с одинарной точностью на GPU, существенно ускоряет работу приложения, т.к. вычисления с двойной точностью на несколько порядков медленнее. Для обеспечения сходимости результатов моделирования диссипативной динамики необходимо наличие большого числа реализаций (M > 103) и генерирования последовательностей случайных чисел непосредственно для каждой точки временной траектории. Формирование массива случайных чисел, необходимого для реализации скачков, выполнялось посредством алгоритма MCG59. 4. Результаты расчетов и полученное ускорение

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

В результате работы над данным вопросом об измерениях было разработано приложение, реализованное на языке CUDA C с использованием пакета CUDA Toolkit 6.5 и технологии MPI. При расчете диссипативной динамики системы с гамильтонианом Вигнера-Вайскопфа (10) использовалась классическая схема Рунге-Кутта четвертого порядка, хорошо подходящая для систем с диссипацией, а в случае происхождения квантовых скачков производилась перенормировка волновой функции по описанным ранее уравнениям ( 13 ).

Для начала были проведены численные эксперименты по вычислению средней энергии нелинейного джозефсоновского осциллятора. Рассмотрим, как ведет себя среднее число фотонов в каждом акте измерения и в среднем по реализациям. Пусть осциллятор в начальный момент времени был заселен на уровень n0 = 28 тогда, если внешнего возбуждения нет, т.е. f0 = 0, среднее число фотонов (средняя энергия) с течением времени уменьшается экспоненциально, что видно на рис. 2, где кривая 1, полученная прямым решением уравнения для матрицы плотности ( 9 ). Однако, данная зависимость описывает лишь поведение системы в среднем, в каждой реализации моменты испускания кванта энергии (фотона) случайны, что было промоделировано на основе квантового метода Монте-Карло (кривые 2, 3 и 4 на рис. 2). При числе реализаций равном M = 100000 наблюдается точное совпадение зависимостей, посчитанных решением уравнения ( 9 ) и квантовым методом Монте-Карло (кривая 1 на рис. 2). При расчетных параметрах ( = 0,005) система должна затухнуть на времени t = 200T = 1256 (ед), что и подтверждается численным расчетом. Кривая 2 для M = 10000 и кривая 3 для M = 1000 квантовых траекторий показывает существенное возрастание погрешности статистического метода, а кривая 4 демонстрирует квантовые скачки в уникальной единичной реализации. 4 1 3 2 Рис. 2. Зависимость среднего числа фотонов n нелинейного осциллятора от времени для разного числа реализаций: кривая 1 для M = 100000, кривая 2 - M = 10000, кривая 3 - M = 1000, кривая 4 - M = 1. Параметры системы: 0  1, f0  0,   0.001,   0.005.

Далее производилось моделирование с возбуждающей силой на интервале времени 0  t  50Т , шаг по времени при расчете скачков dt = 0.02, шаг по времени при расчете методом Рунге-Кутты ddt = 0.000667. Параметры системы для гамильтониана ( 5 ): 0 = 1.02,  = 0.001, f0 = 0.2,  = 0.001.Результаты измерения ускорения работы приложения при 50 уровнях осциллятора, связанных с двухуровневой системой кубитом (т.е. эффективно матрица плотности связанной системы имеет размерность 100х100) при усреднении по различному числу реализаций квантовой системы приведены на рис. 3. Описание параметров запуска приведено в таблице 1. Кол-во GPU 4 Из графика, представленного на рис. 3, видно, что при увеличении числа реализаций (траекторий М в квантовом методе Монте-Карло) наблюдается существенный выйгрыш в ускореении при использовании GPU устройств. Для малого числа реализаций, т.е. если число квантовых траекторий меньше числа потоков на видеокарте, ускорения при использовании одновременно нескольких GPU ускорителей получить не удается (см. на рис. 3 кривые 1, 2, 3). Это связано с тем, что рассылаемое число заданий для вычислений меньше, чем количество потоков даже на одной видеокарте, не говоря уже о нескольких. Для эффективной работы на каждый GPU должно приходиться не менее 1024 путей (и желательно, чтобы количество путей было кратно 1024), т.е. например при M = 4000 достаточно использовать 4 процесса (GPU), ускорения от использования большего количества процессов не будет. Это обусловлено тем, что время расчета одной траектории мало, а при «перегрузке» графического ускорителя наступает эффективная загруженность и достигается оптимальная работа. Прирост производительности вычислений в зависимости от количества используемых графических ускорителей для квантового метода Монте-Карло при условии загруженности (число квантовых траекторий больше, чем число доступных потоков на видеокартах, см. на рис. 3 кривая 4), это говорит о хорошем масштабировании задачи при переходе к множеству GPU. Например, прирост производительности для 100000 квантовых траекторий при использовании 10 GPU увеличился в 6 раз по сравнению с одной видеокартой.

Если сравнивать производительность разработанных параллельных версий программ с использованием технологии CUDA и последовательной версии кода на языке С++ на CPU, то выигрыш в расчетах по времени оказывается очень существенным. Тестирование CPU версии кода производилось на (Intel Core i7 3770K @ 3.50 GHz, использовалось 1 ядро), данные по времени собраны в таблице 2. Видно, что даже для малого числа реализаций (M ~ 100), наблюдается ускорение во времени счета в 5 раз относительно CPU за счет использования графического ускорителя. С ростом числа реализаций GPU версия имеет неоспоримое преимущество, что видно из нижней строки таблицы 2. Таблица 2. Время расчета усредненных M траекторий на GPU устройствах и последовательной версии на CPU. Параметры системы аналогичны тем, которые были взяты для моделирования в таблице 1. Число реализаций, 1 CPU 1 GPU 4 GPU 10 GPU M Тип устройства

M = 1 M = 1000 M = 10000 M = 100000 5. Обсуждение результатов

В работе представлена математическая модель связанной системы "кубит+измерительный" осциллятор. Разработан алгоритм, основанный на стохастическом методе квантовых траекторий (квантовом методе Монте-Карло), для моделирования диссипативной динамики многоуровневой квантовой системы на примере кубита, связанного с квантовым нелинейным осциллятором. Квантовый метод Монте-Карло позволяет осуществить эффективное распараллеливание и программную реализацию на гетерогенном высокопроизводительном вычислительном кластере при расчете процессов релаксации в многоуровневых квантовых системах. Создан программный комплекс, который реализует параллельные вычисления на многопроцессорном кластере (с использованием протокола MPI) и графических процессорных устройств (с использованием технологии CUDA). Продемонстрирована эффективность работы графических ускорителей, масштабируемость и описывается возможность взаимодействия нескольких графических ускорителей с применением технологии MPI. Показано, что применимость GPU ускорителей позволило существенно ускорить вычисления. Развитая в работе техника расчета диссипативной динамики многоуровневых систем естественным образом может быть распространена на более сложные физические системы. Литература 1. Нильсен М., Чанг И., Квантовые вычисления и квантовая информация. М.: Мир, 2006, 824 c. 12. А.В. Боресков, А.А. Харламов, Основы работы с технологией CUDA. Издательство "ДМК Пресс", 2010. 234 c. 15. Гергель В. П. Теория и практика параллельных вычислений. Бином. Лаборатория знаний, 2007. 424 с. Supercomputer simulations method of the coupled qubit and oscillator systems Marina Denisenko, Arkady Satanin, Alexey Liniov, Victor Gergel and Dmitry Pashin Keywords: qubit, non-destructive measurements, quantum trajectory technique(Monte Carlo wave-function method), GPU, CUDA The algorithm based on quantum trajectory technique (Monte Carlo wave-function method) for modeling dissipative dynamics of multi-level quantum system was developed. The simulation of the interaction process a qubit with a quantum nonlinear oscillator based on the developed method has been conducted. It is shown that the measuring result for the average energy (number of photons) of the oscillator gives a possibility to carry out separation (tomography) of qubit states both as an average over an ensemble and as for a single trajectory. Monte Carlo wave-function method allows for effective parallelization on heterogeneous high performance computing clusters in calculation of experimentally observed values in multi-level quantum systems with relaxation processes. The basic principles of parallelism algorithm for solving stochastic Schrödinger equation, which were imple-mented by us in program complex, allowing to carry out high performance parallel calculations using GPU have been described. The efficiency of graphics accelerators, scalability of the method was demonstrated and possibility of interconnection of multiple graphics cards using MPI-technology was described.

References

2. Ladd T. D. , Jelezko F. , Laflamme R. , Nakamura Y. ,. Monroe C & O'Brien J. L. , Quantum computers // Nature, 2010 . V 464, P. 45 .

3. Zagoskin A. M. , Quantum Engineering - Theory Design of Quantum Coherent Structures, (Cambridge: Cambridge University Press), 2011 , 346 p.

4. You , J.Q. , Nori , F. : Atomic physics and quantum optics using superconducting circuits . Nature , 2011 . V. 474, P. 589 .

5. Siddiqi I , Vijay R. , Pierre F. , Wilson C. M. , Metcalfe M. , Rigetti C. , Frunzio L. , and Devoret M. H. , RF-Driven Josephson Bifurcation Amplifier for Quantum Measurement//Phys. Rev. Lett., 2004 . V. 93. P. 207002 .

6. Siddiqi I. et. al., An RF -Driven Josephson Bifurcation Measurements//Phys. Rev. B , 2006 . V. 73. P. 054510 .

8. Plenio M. B. and Knight P. L.The Quantum Jump Approach to Dissipative Dynamics in Quantum Optics// Rev. Mod. Phys., 1998 . V. 70. P. 101 .

9. Vijay R. , Slichter D. H. , Siddiqi I. , Observation of Quantum Jumps in a Superconducting Artificial Atom// Phys. Rev. Lett. , 2011 . V. 106. P. 110502 .

11. Сандерс Д. , Кэндрот Э. , Технология CUDA в примерах . Введение в програмирование гра- фических процессов: Пер. с англ . Слинкина А.А., научный редактор Боресков А.В. М.: ДМК Пресс , 2011 . 232 с.

13. URL: http://hpc-education.unn.ru/files/seminars/_lobachevsky/01_Lobachevsky_general_info. pdf (дата обращения 27.07 . 2015 ).