Русский
Русский
English
Статистика
Реклама

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

Перевод Суперкомпьютеры и клеточные мембраны 2

28.04.2021 10:12:44 | Автор: admin

Предыдущая часть

С самодельным параллельным суперкомпьютером в рюкзаке Клаус Шультен терпеливо ждал в чикагском аэропорту О'Хара, надеясь, что после прибытия из Германии ему не составит труда пройти таможню. Это было летом 1988 года, и Шультен собирался начать новую работу в Университете Иллинойса. В разгар холодной войны, когда напряженность между США и Советским Союзом достигла наивысшего пика, суперкомпьютеры вызывали у администрации Рейгана большой ужас. Хотя Рейган, находясь на своем посту, усилил гонку вооружений и все сопутствующие ей технологические достижения, он хотел, чтобы бурно развивающиеся разработки суперкомпьютеров не попали в руки Советов, которые могли бы создать более совершенное оружие.

Оглавление

  1. Пересматривая устоявшиеся подходы

  2. Предтечи молекулярной динамики

  3. Шультен-числодробитель

  4. Рискованный план

  5. Паяльник вместо экзамена

  6. Транспортировка суперкомпьютера во время холодной войны

  7. Невзгоды теоретической биологии

  8. В поисках сотрудников

  9. Бунт аспирантов

  10. Отказ от платоновской информатики

  11. Ингредиенты NAMD

  12. Награда за декаду

  13. Вычислительный микроскоп достигает совершеннолетия

  14. NAMD в двадцать первом веке


В 1986 году правительство объявило о предполагаемой политике ограничения использования советскими учеными суперкомпьютеров, размещенных в университетах США. Все предельно прозрачно: никакого доступа для ученых из 20 стран коммунистического блока. Но реализация была не столь ясна. Университеты были возмущены этой политикой, беспокоясь о посягательствах на их академическую свободу и исследовательские инициативы, не понимая, как им стать исполнителями политики государственной безопасности.

Фотография самодельного суперкомпьютера, который Клаус Шультен перевозил в рюкзакеФотография самодельного суперкомпьютера, который Клаус Шультен перевозил в рюкзаке

Администрация Рейгана не только собиралась запретить советским ученым проводить моделирование на суперкомпьютерах, но и не хотела, чтобы советские ученые вообще находились рядом, опасаясь, что они могут научиться создавать свои собственные. Но в 1987 году два молодых студента-физика из Мюнхена взяли на себя именно такую миссию по созданию собственного параллельного суперкомпьютера, хотя ни один из них не имел формальной подготовки в этой области. Они не только воплотили свое стремление, но и стоимость их проекта составила около 60 000 долларов, что намного меньше, чем розничная цена Cray-2, популярного суперкомпьютера в конце 1980-х. Их научный руководитель, Клаус Шультен, решил рискнуть всеми грантовыми деньгами, хотя у него не было никаких гарантий, что проект ждет успех, и даже несмотря на то, что он не был экспертом в параллельных вычислениях.

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

Пересматривая устоявшиеся подходы

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

На курсах химии в колледже иногда, помимо учебников, требуются наборы для построения молекул. Студенты собирают молекулы из разноцветных атомов-шариков скрепляя их тонкими палочками. Результатом будет жесткая 3D-модель, которую можно повертеть в руках и рассмотреть со всех сторон. До 1977 года этот статический взгляд на биомолекулу все еще был распространен в некоторых научных кругах. Но постепенно начал происходить сдвиг в концептуальном понимании исследователей, после того, как 1977 год ознаменовался введением вычислительного метода, который прояснил движение атомов, составляющих биологические макромолекулы, что позволяло рассматривать их динамику. Этот сдвиг в сторону использования вычислительных методов для отражения такой динамики стал первым шагом в легитимации вычислительного микроскопа и дал стимул, побудивший Клауса Шультена создать свой собственный суперкомпьютер.

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

В 1968 году Мартин Карплюс был очень успешным химиком-теоретиком, работавшим в Гарварде, и он неожиданно понял, что пришло время вернуться к своей первой любви биологии.

Когда Карплюс был еще мальчиком, отец подарил ему микроскоп, что помогло воспитать в нем любовь к природе. Это пристрастие к природе только усилилось, когда Карплюс посетил лекцию в Бостонской публичной библиотеке под названием "Птицы и их идентификация в полевых условиях", что привлекло его к орнитологии. После, в средней школе, Карплюс был одним из победителей Westinghouse Science Talent Search со своим проектом по семейству чистиковых. Его любовь к биологии продолжалась и в студенческие годы в Гарварде, поскольку в своей автобиографической статье 2006 года он писал: "Я пришел к выводу, что для того, чтобы приблизиться к биологии на фундаментальном уровне (чтобы понять жизнь), необходимо иметь солидный опыт в химии, физике и математике, и поэтому я записался на программу по химии и физике." Однако его аспирантура по биологии не удалась, он перешел на химию и получил докторскую степень осенью 1953 года в Калтехе. (подробности можно найти в автобиографии "Шпинат на потолке: химик-теоретик возвращается в биологию")

Пятнадцать лет спустя, после больших успехов в теоретических исследованиях химии, Карплюс был готов возобновить свой интерес к биологии. Он остановился на Институте Вейцмана в Реховоте, Израиль, как на месте, где он мог проводить время в знаменитой библиотеке и обсуждать идеи в группе Шнейора Лифсона, и взял отпуск из Гарварда, чтобы перестроить свои исследовательские интересы.

Пребывание в Институте Вейцмана в 1969 году дало Карплюсу то, что он искал, и он вернулся в Гарвард с темами для изучения, а именно: происхождение связности гемоглобина, изучение ретинали в зрении и укладка белка.

К началу 1970-х годов один из его учеников, Брюс Гелин, разработал компьютерную программу для вычисления силовых взаимодействий белка просто из его аминокислотной последовательности, а также из координатной разметки его структуры (полученной из экспериментальных данных рентгеновской кристаллографии). Программа была использована в основной работе Гелина по гемоглобину, а затем для изучения ингибитора трипсина поджелудочной железы крупного рогатого скота.

С программой Брюса Гелина к середине 1970-х годов группа Карплюса подготовила почву для применения молекулярной динамики к биомолекулам. К моменту когда к группе присоединился Эндрю Маккаммон, код Гелина обрел законченную форму, реализуя уравнения Ньютона для расчета динамики или отдельных движений атомов белка.

Однако остальные коллеги не очень-то одобряли применение молекулярной динамики к биомолекулам. Как писал Карплюс в своей автобиографии 2006 года: "Когда я обсуждал свои планы с коллегами-химиками, они думали, что такие расчеты невозможны, учитывая сложность точного анализа нескольких атомных систем; коллеги-биологи считали, что даже если бы мы могли сделать такие расчеты, они были бы пустой тратой времени."

Но Карплюса это не разубедило. Он был воодушевлен двумя ответвлениями исследований более простых систем, которые сделали молекулярную динамику возможной для биомолекул. Первой была работа по расчету траекторий нескольких частиц для простых химических реакций, которую Карплюс изучал сам, когда занимался теоретической химией. Эти траекторные расчеты уходили корнями в 1930-е годы. Второе ответвление заключалось в изучении термодинамических свойств большого числа частиц - работа, начатая Берни Олдером и Томом Уэйнрайтом в конце 1950-х годов, когда они впервые применили молекулярную динамику к большой системе взаимодействующих частиц.

Карплюс знал не только о работах 1950-х годов по молекулярной динамике, но и о том, что в середине 1960-х и середине 1970-х годов другие исследователи использовали молекулярную динамику для моделирования жидкостей. Таким образом, с оглядкой на работы предшественников, Карплюс решил рискнуть сделать расчеты молекулярной динамики на биомолекуле, даже если его коллеги по биологии и химии не разделяли энтузиазма. Но преобладающее негативное отношение было не единственным, с чем столкнулись Карплюс и его сотрудники. Расчеты молекулярной динамики, которые они хотели провести, требовали серьезной вычислительной мощности. Карплюс утверждает, что в конце 1970-х годов было трудно получить компьютерное время для выполнения таких больших вычислений в Соединенных Штатах, поскольку большинство больших компьютеров предназначались исключительно для оборонных работ. Но когда возможность принять участие в семинаре в Европе пообещала доступ к компьютеру, Брюс Гелин и Эндрю Маккаммон принялись неустанно работать над подготовкой программы для запуска молекулярной динамики на небольшом белке.

Работа в CECAM (Center Europen Calcul Atomique et Molculaire) на большом компьютере прошла успешно, и троица, Маккаммон, Гелин и Карплюс, вскоре опубликовала свою статью по молекулярной динамике ингибитора трипсина поджелудочной железы крупного рогатого скота. Поскольку молекулярная динамика буквально показывает внутренние движения белка во времени, они смогли воспроизвести 9.2 пикосекунды времени жизни белка. В принципе, программа молекулярной динамики должна была бы воспроизвести время необходимое для развития биологического процесса. И хотя несколько пикосекунд могут показаться небольшим количеством времени для выяснения поведения белка, это послужит одной из направляющих сил в продвижении будущих исследователей к поиску более крупных и быстрых компьютеров. Ученые хотели продлить время моделирования расчета молекулярной динамики, то есть использовать вычислительный микроскоп для длительностей, выходящих за пределы пикосекундного диапазона.

Шультен-числодробитель

Шультен знал о работе Маккаммона, Гелина и Карплюса, потому что он заканчивал аспирантуру в Гарварде, где Карплюс был одним из его преподавателей, а Маккаммон и Гелин были сокурсниками, которых часто можно было встретить в коридоре. Статья 1977 года оказала определенное влияние на Шультена, который был обучен использовать теоретические и математические методы, такие как те, которые используются в химии и физике. "Я понял, что этот вычислительный подход, открывает новые двери для описания проблем, которые нельзя было сделать с чисто теоретическим подходом, который я использовал прежде."

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

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

Несмотря на многочисленные трудности, с которыми пришлось столкнуться Шультену, он был полон решимости использовать молекулярную динамику в своей работе; он был уверен, что это приведет его к новым открытиям, которые будут ценны для науки. "Моя любовь к научным открытиям, заставила меня заняться грязными вычислениями."

Рискованный план

Шультен хотел досконально разобраться в фотосинтезе, при котором солнечный свет превращается в биохимически пригодную энергию. Но важной макромолекулой в этом процессе является большой белок, который находится внутри мембраны. Имитировать белок сам по себе и пренебрегать мембраной и жидкостью, которая его окружает, было не желательно, потому что такая изоляция не является естественной средой для белка. Этот белок состоит примерно из 12 000 атомов, а мембрана и вода, окружающие его в естественной среде, добавляют еще 100 000. В конце 1980-х годов ни один суперкомпьютер даже близко не был способен справиться с этой задачей. Таким образом, Шультен решил сосредоточиться на понимании и моделировании только части мембраны в воде, когда выпала возможность выполнить расчет на Cray-XMP, но это позволило воспроизвести процессы всего в несколько пикосекунд, что натолкнуло на мысль: нужен свой суперкомпьютер на котором можно гонять расчеты год или больше, чтобы действительно понять механизм.

В то время как Шультен, преподававший физику в Мюнхенском техническом университете, задавался вопросом, как завести себе суперкомпьютер, молодой студент того же университета задавался вопросом, как получить возможность быстрый компьютер построить. На самом деле, теоретически, Гельмут Грубмюллер знал, как сделать компьютер быстрее, но ему было интересно, сможет ли он заставить кого-то другого оплатить это. При поступлении в Мюнхенский университет в 1985 году, двадцатилетний Грубмюллер осознал, как выбор специальности определил его цели: "Когда я начал заниматься физикой, я быстро понял, что для того, чтобы узнать больше о том, как работает природа, может потребоваться мощная вычислительная машина. В конечном счете, именно поэтому я был так очарован компьютерами; вы могли делать вычисления, которые иначе не осилить."

Во втором или третьем семестре он спаял многопроцессорный компьютер, который будет намного быстрее, чем то, что было доступно обычным студентам в те дни. Хотя технически это был параллельный компьютер, Грубмюллер просто назвал его мультипроцессором. Это была середина 1980-х, и не было никакой Всемирной паутины, чтобы проконсультироваться по техническим вопросам. Вместо этого он читал любые книги о микропроцессорах, какие только мог найти, и разговаривал по телефону с представителями компаний о технических характеристиках продаваемых ими деталей. В качестве наиболее важных для продвижения вперед источников информации для него были технические спецификации чипов, например, из линейки Motorola 68 000, которые он использовал в качестве процессора. Это увлечение быстро опустошило бюджет студента.

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

Клаус Шультен в то время преподавал потоку Грубмюллера курс механики, и он помнит тот день, когда молодой человек подошел к нему со своим самодельным многопроцессорным компьютером. "Он открыл свою сумку, достал оттуда суперкомпьютер, показал мне результаты работы и рассказал о характеристиках, и я чуть не упал со стула", Шультен был впечатлен демонстрацией некоторых возможностей мультипроцессора. "Поэтому я сразу же включил его в свою группу, сказав ему, что в группе есть и другие люди, которые тоже много знают о компьютерах."

Одним из таких людей был Гельмут Геллер, который, как впоследствии узнал Шультен, пользовался некоторой славой в университете. Еще в старших классах Геллер начал сам изучать компьютеры. Он был счастливым обладателем одного из первых персональных компьютеров, PET 2001, представленный Commodore в конце 1970-х. Там было всего 8 килобайт оперативной памяти и кассета для архивирования данных. "Я начал с изучения BASIC с отладкой программ на этой машине. Потом же, начал писать машинный ассемблерный код", вспоминает Геллер. У него в родном городе было несколько друзей, которые тоже интересовались компьютерами и учили друг друга.

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

Клаус Шультен, Гельмут Грубмюллер, Гельмут Геллер и одна из плат самодельного суперкомпьютера,1988.Клаус Шультен, Гельмут Грубмюллер, Гельмут Геллер и одна из плат самодельного суперкомпьютера,1988.

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

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

"Иногда приходится рисковать", оправдывает Шультен свое решение позволить Грубмюллеру и Геллеру построить параллельный суперкомпьютер. Риск для Шультена был огромен. Хотя у него был бы суперкомпьютер для моделирования молекулярной динамики, если бы это каким-то образом не сработало, ему было бы нечего показать. В то время он не был экспертом в области параллельных вычислений и имел мало сведений для оценки осуществимости проекта. Кроме того, ему придется объяснять финансирующему агентству, как он вложил десятки тысяч немецких марок в самодельный компьютер, тогда как деньги предназначались для оборудования, которое обычно поставлялось с гарантиями и контрактами на обслуживание.

Паяльник вместо экзамена

Геллер и Грубмюллер, однако, составили отличную команду, и то, что они пытались сделать, было немалым подвигом. Им нужно было не только построить суперкомпьютер, но и написать программное обеспечение для параллельной работы на нем. "У нас с Гельмутом Грубмюллером были очень похожие, но в то же время комплиментарные компьютерные знания. Он больше занимался аппаратным обеспечением, а я программным.", вспоминает Геллер. Однако они разрабатывали и аппаратное, и программное обеспечение вместе, консультируясь друг с другом на каждом шагу.

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

В конечном итоге компьютер был собран и назван T60. За основу использовался кластер процессоров, известных как транспьютеры, потому что якобы из них можно было так же легко собрать схему, как из блоков Lego. В результате вышло 60 узлов: десять самодельных плат, по шесть транспьютеров на каждую. Все компоненты были спроектированы, спаяны и собраны вместе: транспьютеры, оперативная память, источник питания, вентиляторы (аэраторы), корпус, шины питания, печатные платы и крепления.

Поскольку пайка каждой из десяти печатных плат занимала много времени, к работе подключились другие члены исследовательской группы Шультена. Каждый участник получал детали и лист инструкций и начинал паять плату, что занимало по крайней мере около шести недель, а затем подписывал на ней свое имя. "Это действительно показывало исследовательскую личность студента", утверждает Шультен. "Работа, которую каждый в группе должен был проделать, действительно очень близко отражала то, как они на самом деле будут вести себя в команде и позже в науке." На самом деле Шультен считает, что вся эта пайка была даже лучше, чем традиционные результаты экзаменов, в определении того, какие кандидаты будут хорошими аспирантами, чтобы принять их в группу.

Геллер вспоминает о трудностях, с которыми столкнулась команда при написании кода для параллельной машины, кода, который они назвали EGO. "Я помню, что когда я впервые сделал параллельные вычисления, которые хорошо масштабировались, я был очень взволнован, но вскоре я обнаружил, что все результаты были неправильными! У меня была быстрая параллельная программа, дающая неправильные результаты." Найти причину этого было нелегко. Шультен купил кусок промышленного марципана, до которого Геллер был большой любитель, и время от времени отрезал ему ломтики, чтобы побудить его продолжать. После долгих добродушных насмешек со стороны других членов группы за его быстрые, но неправильные вычисления, Геллер в конце концов понял, что ему нужно правильно синхронизировать параллельные задачи, так как сначала он дал им слишком много свободы; решение заняло много дней и долгих вечеров работы.

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

Проект по созданию и программированию Т60 начался в конце 1987 года. На пол пути к его завершению Шультен получил новую работу в Соединенных Штатах, в Университете штата Иллинойс в Урбана-Шампейн. Вместо того чтобы отправить Т60 в Соединенные Штаты и неделями ждать, пока он прибудет, а затем пройдет таможню, Шультен решил просто пронести компьютер в рюкзаке.

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

Существовал специальный список ограничений и запрещенных материалов для экспорта в коммунистические страны, который был составлен под эгидой COCOM, комитета, созданного после Второй мировой войны для ограничения поставок оружия в Советский Союз. Список был единогласно согласован Соединенными Штатами, союзниками по НАТО и Японией. В 1984 году, после долгих ожесточенных споров между Соединенными Штатами и Западной Европой, особенно по поводу экспорта компьютеров, COCOM пересмотрела свой список товаров и в основном разрешила свободно экспортировать персональные, но не суперкомпьютеры.

Администрация Рейгана не только хотела ограничить бизнес, как внутренний, так и зарубежный, от экспорта суперкомпьютерных технологий в СССР, она также хотела держать советских ученых подальше от своих недавно созданных суперкомпьютерных центров. В 1985 году Национальный научный фонд выделил четырем университетам 200 миллионов долларов на открытие суперкомпьютерных центров в их кампусах. Этот шаг был признан необходимым и санкционирован Конгрессом, чтобы предоставить академическим исследователям доступ к суперкомпьютерам, которые обычно ограничивались оборонными работами, выполняемыми Пентагоном и Национальными лабораториями. Как сказал один из ключевых ученых, лоббировавших правительство США, использование стандартных компьютеров вместо упрятанных за бюрократическими препонами суперкомпьютеров "похоже на езду на лошади и возу, когда над головой летают самолеты". Но в 1986 году администрация Рейгана предложила запретить доступ к этим суперкомпьютерным центрам советским ученым. Таков был климат летом 1988 года, когда Клаус Шультен решил взять свой самодельный суперкомпьютер на трансатлантический рейс и пронести его через таможню в Чикаго.

Шультен знал, что из-за скорости обработки данных компьютер Т60 должен быть зарегистрирован, особенно при экспорте или импорте. Он подошел к таможеннику и показал ему суперкомпьютер.

Что это? спросил таможенник.

Это компьютер, ответил Шультен.

Зачем Вы мне его показываете?

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

Затем Шультен объяснил офицеру некоторые технические детали машины. Наконец таможенник посмотрел на Шультена и сказал: "Пожалуйста, положите его обратно и ступайте."

Невзгоды теоретической биологии

Когда Шультен прибыл в Иллинойс, его ученики закончили собирать и программировать свой Т60. Программа, получившая название "EGO", была написана на языке OCCAM II языке, физическим воплощением которого, по сути, и являются транспьютеры. Точно так же, как создание T60 требовало больших усилий, написание параллельного кода было столь же трудной задачей. Одним из ключевых показателей в параллельном программировании является то, насколько хорошо масштабируется программное обеспечение; хорошее масштабирование (или параллелизуемость) означает, что время выполнения вычислений приближается к ускорению кратному количеству логических потоков. Например, сбор яблок человеком хорошо параллелизуется: по мере того, как все больше людей участвует в сборе, время выполнения работы уменьшается.

Мембранная структура, смоделированная на Т60 в течение двадцати месяцев, охватывала 24 000 атомов. Изображение из Heller et al., 1993.Мембранная структура, смоделированная на Т60 в течение двадцати месяцев, охватывала 24 000 атомов. Изображение из Heller et al., 1993.

Было очень важно, чтобы EGO хорошо масштабировалось, потому что расчет мембраны, который Шультен и его группа хотели запустить, был очень большой системой. На самом деле расчет занял двадцать месяцев безостановочных вычислений. Система, которую они изучали, мембрана из липидов, окруженная водой, состояла из 23 978 атомов. Для сравнения, число атомов в системе биомолекул, которую Маккаммон, Гелин и Карплюс изучали в 1977 году, было на два порядка меньше.

После таких невероятных усилий, приложенных к проблеме мембраны, можно представить себе удивление Шультена, когда их с Гельмутом Геллером и Михаэлем Шефером статья была отклонена в течение недели, хотя результаты хорошо согласовались с экспериментом. Статья содержала результаты теоретической биофизики. Шультен признает, что опубликовать теорию, связанную с биологией, очень сложно, особенно в журналах с высоким импакт-фактором. Он говорит, что даже сегодня у него больше шансов получить признание, если он объединится с экспериментаторами. В своей автобиографической статье 2006 года Мартин Карплюс суммирует препятствия, с которыми сталкиваются теоретики при публикации результатов, связанных с биологией: "Как и раньше, сегодня эта проблема все так же распространена, то есть, если теория согласуется с экспериментом, она не интересна, потому что результат уже известен, тогда как если кто-то делает прогноз, то он не может быть опубликован, потому что нет доказательств того, что прогноз верен."

Тайна отказа так и не была раскрыта. "Я до сих пор не знаю, почему ее отвергли", рассказывает Шультен. "Я был так взбешен, что, хотя я опубликовал много статей в этом журнале, я сказал им, что больше никогда не буду там публиковаться." Геллер, Шефер и Шультен в конце концов опубликовались в Journal of Physical Chemistry, и статья сейчас высоко цитируется.

Этот набег на молекулярную динамику такой гигантской для своего времени системы дал обоснование для использования параллельного компьютера в качестве вычислительного микроскопа. Было важно взять для анализа большой сегмент мембраны, потому что реальные клеточные мембраны непрерывны и не имеют края как такового. Другие типы обычных микроскопов не могли уловить некоторые критические особенности мембраны, которые проясняла молекулярная динамика. "Когда вы смотрите на мембрану в течении какого-то времени, то замечаете, что она вся текучая и движущаяся. Что довольно трудно рассмотреть с помощью микроскопов. А мы действительно видели с помощью компьютера сам процесс и могли провести сравнение с большим количеством экспериментов. Они измеряют некие усредненные характеристики, которые придавали нам убеждение в том, что то, что мы действительно видели с помощью компьютера, было правильной картиной." Успех вычислений только подогрел аппетит Шультена к изучению все больших и больших систем. Шультен не понимал, что его жажда массивных систем на параллельных машинах приведет к своего рода студенческому бунту и программному продукту под названием NAMD, который определит его карьеру.

В поисках сотрудников

В 1989 году Шультен основал группу теоретической биофизики [Theoretical Biophysics Group] (которая впоследствии стала группой теоретической и вычислительной биофизики) в Институте Бекмана Иллинойского университета. В 1990 году он получил двухлетний грант от Национального института здравоохранения, которого более или менее должно было хватить, чтобы центр заработал. К этому времени Шультен лучше разбирался в параллельных вычислениях, по сравнению с теми временами, когда он вложил все свои деньги в параллельную машину еще в Мюнхене; он понял, что может изучать массивные биомолекулы, используя преимущества параллельных машин. Уже существовали программные коды, некоторые даже в свободном доступе, которые занимались молекулярной динамикой, но Шультен понял, что его потребности намного превышают их возможности.

На самом деле, Шультен понял, что его потребности превышают даже его собственные возможности как вычислительного биофизика, если он хочет осуществить свой план по изучению более массивных биомолекул в их естественной среде. Шультен решил выйти за пределы своего поля. "На мой взгляд, это было некоторое предвидение с его стороны, начать искать опытных информатиков", вспоминает Лаксмикант Санджай Кале, один из ученых, к которым обратился Шультен. Мотивированное тремя причинами, суждение Шультена вынудило его обратиться за помощью к компьютерщикам в своем университете.

Во-первых, до 1991 года программированием молекулярной динамики в группе Шультена всегда занимались его студенты-физики. После того, как оба Гельмута написали EGO для самосборного компьютера, другой студент, Андреас Виндемут, написал код по молекулярной динамике для Connection Machine.

Как уже упоминалось ранее, в 1985 году в Соединенных Штатах были созданы четыре центра для предоставления академическим исследователям доступа к суперкомпьютерам. Четыре центра находились в Корнелле, Принстоне, Калифорнийском университете, Сан-Диего и Иллинойском университете, в Шампейн-Урбане, где находился Шультен. Иллинойский центр назывался NCSA, и группа Шультена имела доступ к Connection Machine, которая размещалась там. Видя программы для T60 и для Connection Machine Шультен понял, что ему нужны навыки разработки программного обеспечения в группе, чтобы убедиться, что программное обеспечение написано таким образом, чтобы в будущем оно было понятно и открыто для модификаций. Его студенты-физики просто не были обучены таким навыкам.

Во-вторых, Шультен понял, что методы программирования меняются, и особенно то, что FORTRAN, разработанный в 1950-х годах, возможно, не будет лучшим вариантом, поскольку начали приобретать известность другие языки. "Нам нужны были лучшие языки программирования, которые были бы более систематичными".

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

Имея в виду вышеупомянутые причины и необходимость получить более солидный грант на обновление своего центра, грант, который обеспечит финансирование на пять лет, а не только на два года, Шультен обратился к двум ученым-компьютерщикам из университета Иллинойса и попросил их принять участие в работе группы. Один из них, Боб Скил, был экспертом в области численных алгоритмов, а другой, Санджай Кале, был экспертом в параллельном программировании. В 1992 году Шультен получил пятилетний грант от Национального института здравоохранения и таким образом приступил к основам программного кода под названием NAMD.

Бунт аспирантов

На самом деле Шультен воспринимает раннюю работу над T60 как приготовление почвы для NAMD это вкупе со студенческим бунтом, который произошел в начале 1990-х годов. В это время Шультен надеялся продолжить использовать коды молекулярной динамики, которые были разработаны для Connection Machine и Т60. Аспирант Боба Скила, Марк Нельсон, изучал код Connection Machine, а аспиранты Шультена, Билл Хамфри и Эндрю Далк, изучали код Т60. Все трое были непомерно раздражены своей работой. "Я бился головой об этот код в течение нескольких месяцев, рассказывает Нельсон, и никуда не продвигался. Комментариев не было, и все имена переменных были сокращениями немецких названий, что как-то мне не помогало."

Трое аспирантов, Нельсон, Хамфри и Далк, решили, что их задача будет проще, если они просто напишут новый код с нуля. Они тайно работали около месяца, пока не получили жизнеспособный код. Взяв на себя большой риск перед своими научными руководителями, они представили на собрании группы идею создать совершенно новый код и оставить в покое предыдущие реализации, которые было бы сложно поддерживать. "Они хотели написать все на C++ и боялись, что профессор пристрелит их за дерзкий отказ от работы с проверенным кодом." смеется Шультен. Однако на самом деле, он был в восторге от этой инициативы. "Так что я сам был очень взволнован, потому что увидел здесь ситуацию, которая была почти такой же, как у Гельмутов до них. Студенты осознавали риск и были готовы взять на себя ответственность."

На тот момент группа Шультена уже была готова создать программный код с нуля, потому как в команде появились настоящие компьютерщики. Один из них, Санджай Кале, попал в проект в следствие случайного стечения обстоятельств и вскоре узнал, что ему будет нужна толстая кожа для того, что его ждет.

Отказ от платоновской информатики

В начале 1980-х годов, будучи аспирантом в SUNY Stony Brook в Нью-Йорке, Кале начал работать над параллельным логическим программированием. В частности, он использовал язык Пролог для изучения искусственного интеллекта. "Так уж получилось, вспоминает Кале, что если вы осмотритесь в этой области сейчас, то увидите огромное количество людей, которые в середине 1980-х годов кодили на Прологе."

Санджай Кале, специалист по компьютерам в Университете Иллинойса с 1985 года.Санджай Кале, специалист по компьютерам в Университете Иллинойса с 1985 года.

Когда он получил работу в Университете Иллинойса, он продолжал работать над параллельным логическим программированием и получил должность в 1991 году. Пока он занимался этим более или менее "чистым" исследованием, основанным на компьютерных науках, он начал обдумывать применение своей работы. "Я занимался поиском на основе состояний, искусственным интеллектом в общем, и его распараллеливанием в частности, задачами а ля коммивояжера и т. д.", рассказывает Кале об этом периоде времени. "Но мои интересы начали смещаться в сторону: как мы можем применить это к задачам, которые есть у инженеров и ученых?"

По удивительному совпадению, сразу после того, как Кале получил должность и собирался отправиться в академический отпуск в Индию на семестр, к нему подошел Клаус Шультен, чтобы спросить, будет ли Кале участвовать в гранте 1991 года в качестве эксперта по параллельным вычислениям. У Кейла внезапно появился прикладной проект, буквально на его собственном заднем дворе. Он провел часть своего отпуска, читая о биологии и биологическом моделировании, а когда вернулся из Индии, начал изучать две программы, которые Шультен использовал для параллельных вычислений молекулярной динамики.

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

Во-первых, он понял, что компьютерщики часто занимаются абстрактными исследованиями. "Это не работает, потому что вы не концентрируетесь на реальных проблемах", говорит Кале. Поэтому я называю это платоновскими компьютерными науками. Идея выглядит красиво, поэтому вы ее развиваете, а затем перебрасываете ее через забор в реальный мир. А реальному миру все равно, потому что он занят своими проблемами."

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

К началу 1990-х годов Кале работал над прикладными проектами: Клауса Шультена по молекулярной динамике и еще одним по гидродинамике с сотрудником в области машиностроения. Что касается молекулярной динамики, то Кале изучал предыдущие коды, разработанные группой Шультена. Примерно в 1994 году аспиранты предложили начать все сначала. "Мы решили явно оформить его как программу", вспоминает Кале. "Сделать шаг назад и разработать его как программу, а не позволять ему быть органически выращенным." Этот формальный дизайн привел к прототипу NAMD. Название первоначально было акронимом "Not Another Molecular Dynamics", который придумал аспирант Билл Хампри, вдохновленный названием компилятора аналогичной природы; вскоре аббревиатура стала обозначать "NAnoscale Molecular Dynamics". Явное требование NAMD быть параллельным с самого начала будет иметь долгосрочные последствия с точки зрения успеха программного обеспечения.

Хотя Кале исполнил свое желание работать над приложениями, последнее десятилетие века сопутствовалось тяжелым трудом. Кале и его аспиранты разработали язык параллельного программирования, названный ими Charm++, который стал очень важным и ценным для успеха NAMD. Его собственные исследования в области компьютерных наук в течение десятилетия были сосредоточены на этом. "То, что мы предлагали, было параллельным C++ или параллельным объектно-ориентированным языком, который был по-своему новым и существенно отличался от того, как люди занимались параллельным программированием", рассказывает Кале. "Так что это был мой вызов заставить людей принять это."

Этому предприятию также характерно большое число поражений. У Кале остались стопки отклоненных грантовых предложений в качестве напоминания. Рецензенты говорили, что то, что он предлагал, либо невозможно, либо устарело. Его публикации тоже не всегда ценились: "Мы писали статьи, и я говорил: "Это продемонстрировано в молекулярной динамике". А люди отвечали: "Это просто параллельный C++. Многие люди делают параллельный C++". Они не обращали внимание на элемент новизны. И работы продолжали игнорироваться".

Кале также считает, что работа на чужом поле была еще одной причиной, по которой ему приходилось бороться. "Люди говорят о междисциплинарных исследованиях, но на самом деле это очень трудно реализовать. Они не очень уважают междисциплинарные исследования ни в одной из этих дисциплин." Оглядываясь назад на то время, Кале, как и Шультен, оценивает траекторию своей профессиональной жизни после того, как в 1990-х годах он принял решение перейти к практическим задачам: "Я все еще не советую молодым людям делать то, что сделал я, если они действительно не понимают, какой риск таит этот карьерный путь." Пройдут годы, прежде чем этот риск окупится для Кале.

Ингредиенты NAMD

Многие факторы формировали направления в которых развивался NAMD в 1994 году, не последним из которых было то, что теперь он хорошо финансировался грантом NIH. Аспиранты настаивали на использовании языка С++, относительного новичка в области программирования. "C++ предложил хорошую комбинацию объектно-ориентированного языка, заявляет бывший студент Марк Нельсон, но мы все равно могли бы написать ключевые вычислительные части на C и получить необходимую скорость." Использование C++ в NAMD сделало его уникальным среди солверов молекулярной динамики того времени.

Поскольку NAMD должен был быть адаптирован для параллельных машин, разработчики должны были учитывать, как многие процессоры, составляющие параллельную машину, будут взаимодействовать друг с другом при выполнении вычислений молекулярной динамики. Кале и Аттила Гурсой, его аспирант, работали над средой Charm++ и кое-чем, что называлось message-driven execution. Хотя технические детали выходят за рамки этой статьи, достаточно сказать, что Charm++ использовал message-driven execution, чтобы эффективно позволить процессорам общаться друг с другом во время выполнения молекулярной динамики.

В то время как C++ и Charm++ были уникальными элементами для NAMD, параллельный дизайн с нуля был еще одной отличительной чертой. В начале 1990-х годов другие очень успешные коды молекулярной динамики, такие как Amber и CHARMM (который является гарвардским продуктом и не следует путать его с Charm++, который был разработан в Иллинойсе), приняли другую тактику, когда стало очевидно, что программное обеспечение должно быть разработано для параллельных компьютеров.

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

Санджай Кале приводит преимущества, которые дизайн с нуля давал NAMD: программное обеспечение хорошо масштабировалось и могло обрабатываться любыми новыми машинами, которые быстро появлялись на рынке. "Этот дизайн выдержал испытание временем и до сих пор актуален", размышляет Кале. "В значительной степени архитектура параллельной программы осталась прежней, хотя в то время мы работали на кластере HP, с восемью процессорами. А теперь мы работаем на 200 000 с лишним процессорах.

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

Награда за декаду

Почти через десять лет после того, как Кале начал свою одиссею в мир параллельных вычислений, он, наконец, был вознагражден. Каждый год проводится международная конференция Supercomputing. Премия Гордона Белла вручается ежегодно на конференции в знак признания достижений в области параллельных вычислений. И в 2002 году она была присуждена разработчикам NAMD.

Хотя Кале, возможно, работал в течение многих лет над такими прикладными проектами, как молекулярное моделирование и гидродинамика, он также занимался более или менее чистыми исследованиями в области компьютерных наук. Он считает, что премия Гордона Белла принесла его группе заслуженную славу. "На самом деле, я очень благодарен за это, отмечает Кале, потому что наши идеи в области информатики остались бы недооцененными. И было бы трудно добиться их признания, если бы не успех NAMD". Из-за своей репутации в NAMD у Кале теперь есть проекты в области моделирования погоды, вычислительной астрономии и вычислительной химии.

Вычислительный микроскоп достигает совершеннолетия

К этому времени вычислительный микроскоп стал более сложным. В 1977 году было проведено первое моделирование молекулярной динамики на небольшом белке в вакууме, которое отображало 9,2 пикосекунды. К 1996 году NAMD моделировал систему с 36 000 атомами, рецептор эстрогена с сегментом ДНК в соленой воде, в течение 50 пикосекунд. Расчет выполнялся в течение двух-трех дней на 8-процессорном кластере. В 2004 году для системы аквапорин-1, состоящей из 81 065 атомов, расчет проводился в течение 22,4 часов на параллельной машине с 128 процессорами и отобразил 5 наносекунд времени жизни системы.

Эти достижения утвердили вычислительный микроскоп для научного сообщества. "В прежние времена компьютерное моделирование не было таким точным", вспоминает Шультен. Поэтому было бы высокомерием просто сказать: 'Это вычислительный микроскоп', люди посмеялись бы над нами и сказали бы: 'У вас очень плохой микроскоп'.

Когда его спрашивают, какие силы заставили вычислительный микроскоп совершенствоваться, Шультен ясно понимает, что есть более глубокая причина, чем очевидная, стоящая за его мотивом постоянно совершенствовать NAMD для создания более крупных систем. "Мы стали больше не потому, что я хотел бить себя в грудь и говорить, что я могу сделать больше атомов, чем остальные, а скорее по ясной интеллектуальной причине. И причина в том, что я думаю, что для самой природы живых систем, в частности живых клеток, необходимо описать, как биологические макромолекулы, в частности белки, собираются и взаимодействуют. Ключевым моментом является просто то, что эти макромолекулы должны найти друг друга в нужной пропорции, в правильной геометрии, в нужном месте в клетке, и они должны сотрудничать." И для того, чтобы описать это, NAMD пришлось усовершенствовать до такой степени, что он мог описать десятки, сотни или тысячи белков.

Часто вычислительный микроскоп правильно предсказывал то, что наблюдал эксперимент; иногда он опережал эксперимент. Одним из таких случаев был пример использования NAMD для запуска моделирования "управляемой молекулярной динамики". В управляемой молекулярной динамике можно применить внешние силы к белку и отобразить его последующее поведение во времени. Анкириновый повтор это определенная последовательность аминокислот, и эта последовательность содержится в более чем 400 белках человека. В 2005 году Маркос Сотомайр, Дэвид П. Кори и Клаус Шультен изучали упругие свойства этих анкириновых повторов с помощью управляемой молекулярной динамики, в конечном счете, чтобы лучше понять слух и чувство равновесия; год спустя, в 2006 году, эксперимент подтвердил их результаты моделирования. В высоко цитируемой научной статье 2007 года Одномолекулярные эксперименты in Vitro и in Silico Сотомайр и Шультен утверждают, что эксперименты in silico "стали мощным инструментом, дополняющим и направляющим эксперименты in vitro". Термин 'in silico' относится к исследованию, выполненному с помощью компьютерного анализа. Обоснование вычислительного микроскопа как мощного инструмента только усиливалось по мере того, как этот эксперимент предвосхитил результаты эксперимента in vitro.

NAMD в двадцать первом веке

NAMD появился на свет в начале 1990-х годов, и семя его восходит к самодельному параллельному компьютеру, построенному в конце 1980-х годов, и программе EGO, которая была написана для этого конкретного суперкомпьютера. Шультен хотел управлять все большими и большими системами, и поэтому NAMD родился, а затем вырос, потому что он хорошо финансировался. Ключевым игроком, стоящим за успехом NAMD, является ведущий разработчик Джим Филлипс, бывший аспирант Шультена в 1990-х годах. Теперь Шультен называет Филлипса "отцом NAMD", поскольку Филлипс разрабатывал творческие решения для параллельного программного обеспечения почти два десятилетия.

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

Само собой разумеется, что исследования вирусов помогут в борьбе с болезнями. В 2009 году Всемирная организация здравоохранения объявила пандемию свиного гриппа H1N1pdm. Было быстро обнаружено, что этот штамм гриппа приобрел устойчивость к Тамифлю, лекарству, которое лечит людей, недавно инфицированных гриппом.

Коллапс пустого капсида. Три снимка сверху взяты из моделирования всех атомов , а шесть внизу-из крупнозернистого моделирования.Коллапс пустого капсида. Три снимка сверху взяты из моделирования всех атомов , а шесть внизу-из крупнозернистого моделирования.

Кроме того, было обнаружено, что вирус птичьего гриппа H5N1 также обладает устойчивостью к Тамифлю. Используя NAMD, а также молекулярную динамику и управляемую молекулярную динамику, сотрудники Университета Иллинойса и Университета Юты объединили усилия, чтобы найти основу для этой лекарственной устойчивости. Их открытие, согласно их публикации 2010 года, "предполагает, как мутации нарушают связывание лекарств и как новые лекарства могут обойти механизмы устойчивости."

Вирус табачной мозаики один из самых маленьких известных вирусов; этот факт в сочетании с достижениями в области компьютерных технологий позволил моделировать эту форму жизни с помощью NAMD целиком. Исследователи из Иллинойса объединились с исследователями из Калифорнийского университета в Ирвине, чтобы изучить стабильность вирусной частицы, а также ее составных частей. Результаты, опубликованные в 2006 году, были первым полностью атомарным моделированием живого организма, охватывающим более миллиона атомов. Полученные результаты установили, какие факторы имеют решающее значение для структурной целостности всего вируса, а также какие факторы могут направлять сборку частицы.

Вирусы это очень примитивные частицы, и многие из них состоят только из белковой оболочки, или капсида, окружающего нить ДНК или РНК. Для размножения вирус проникает в клетку-хозяина и захватывает ее механизмы, чтобы создать больше вирусных частиц. В то время как капсид играет роль защиты вируса, он также должен каким-то образом стать нестабильным и высвободить свои внутренние компоненты в клетку-хозяина для размножения. Распознавание движения белков, составляющих капсид во время такого высвобождения, является приоритетной целью, и исследование, проведенное в 2006 году, изучало динамику различных капсидов. Относительно новый метод крупнозернистая молекулярная динамика, позволял моделировать целые капсиды, которые до этого времени были недоступны. Ученые из Университета Иллинойса усовершенствовали метод крупнозернистой молекулярной динамики и использовали NAMD для иллюстрации стабильности капсида для нескольких вирусов. С помощью новой техники стало возможным моделировать капсиды размером более 10 нанометров и наблюдать переходы между стабильными и нестабильными структурами. Эта работа может дать полезную информацию для борьбы с вирусными заболеваниями.

Будущее для NAMD выглядит радужным. Цель Клауса Шультена понять биологическую организацию, саму суть того, что делает клетку живой, маячит на горизонте после сорока лет самоотверженности и творчества. "Так что теперь я могу потихоньку заняться тем, о чем мечтал, замечает Шультен о биологической организации, ближе к отставке. И другие люди следующие за мной могут использовать наработки сейчас и могут делать это намного лучше. И это, конечно, тоже делает меня счастливым." И вычислительный микроскоп продолжает прояснять все большие и большие системы. "Сегодня у нас есть симуляция на 20 миллионов атомов", хвастается Шультен. "И у нас есть симуляция, чтобы быть готовыми к следующему поколению компьютеров, которые потянут 100 миллионов атомов." Сочетание понимания и принятия риска, безусловно, окупилось с точки зрения научных открытий. От выяснения рецептора эстрогена, чтобы лучше понять рак молочной железы, до освещения вирусов, чтобы лучше бороться с болезнями, вычислительные проекты, ставшие возможными с помощью параллельных компьютеров, будут продолжать раздвигать границы. И до тех пор, пока смелые ученые идут на большой личный риск, пытаясь сделать то, что никогда не делалось раньше, и продолжают переосмысливать новые пути для того, чтобы компьютер действовал как научный инструмент, помогающий исследователям, будущее вычислительного микроскопа, кажется, наверняка принесет больше открытий.

Продолжение следует...

Подробнее..

Перевод Суперкомпьютеры и клеточные мембраны 3

04.05.2021 12:13:46 | Автор: admin

Предыдущая часть

Оглавление

  1. Они не станут обманывать невинную молодежь

  2. Аквапорины, давно искомые водные каналы

  3. Пращи и стрелы

  4. Вход в новый класс белков: механочувствительные каналы

  5. Очарование воздушного шара

  6. Методологическая разработка для Энигмы

Они не станут обманывать невинную молодежь

В 1992 году Клаус Шультен столкнулся с большой дилеммой, за которую был бы благодарен любой ученый. Как лучше потратить полмиллиона долларов на оборудование? Дело в том, что профессор только что узнал, что ему уготовано финансирование от Национального Института Здравоохранения (NIH) для создания центра вычислительной биологии. Сначала он выдержал двухлетний испытательный срок, и когда все прошло хорошо, в 1992 году ему дали стартовые средства на пятилетний период. Шультен считает, что по большей части успех был связан с самодельным суперкомпьютером, созданным его группой, который ясно продемонстрировал роль вычислений в биологии.

В 1985 году Национальный научный фонд санкционировал запуск четырех центров с суперкомпьютерами по всей территории Соединенных Штатов, что должно было дать рядовым исследователям доступ к этим мощным машинам. Один из них оказался в Иллинойском университете, где работал Шультен, когда в 1992 году получил зеленый свет на создание своего центра финансируемого NIH. На самом деле, наличие суперкомпьютерного центра Урбана-Шампейн является одной из причин, по которой Шультен устроился на новую работу в США. Ему нравился акцент на вычислениях в Иллинойсе.

Когда он рассказал о своей проблеме тамошним суперкомпьютерщикам, то сразу получил от них предложение. Используйте свои деньги, чтобы помочь нам купить суперсовременный суперкомпьютер под названием Connection Machine; у вас будет доступ к нему, но нам также понадобятся две из ваших пяти позиций в центре NIH, чтобы облегчить сделку. Шультен никогда не смог бы позволить себе купить собственную Connection Machine, поскольку она стоила более пяти миллионов долларов.

Шультен рассматривал это как вариант. Его группа уже имела опыт использования Connection Machine, чтобы провести кое-какие исследования нейронных сетей. Но ему нужно было потратить свои деньги на машину для моделирования молекулярной динамики совершенно отличной от исследований нейронных сетей. Шультен осмыслял все, что он знал о Connection Machine из разговоров со многими людьми. Он всерьез рассматривал Connection Machine как выбор, но был настроен скептически. "У меня были большие подозрения в отношении этих компьютеров, потому что они всегда описывались с очень эзотерической точки зрения, а не с точки зрения практики, то есть особо не рассказывалось о том, что вы действительно можете вычислить с их помощью."

Оглядываясь назад, Шультен резюмирует дилемму, с которой он столкнулся в 1992 году: "Я хочу получить лучший компьютер для вычислительной биофизики, но какой для этого подойдет, должен разобраться сам." Поэтому он инициировал разведывательную миссию. Он отправил горстку студентов по всей стране в компьютерные компании и институты чтобы поспрашивать, в какой компьютер лучше инвестировать. Сам Шультен специально туда не ходил он чувствовал, что компьютерные компании увидят в его учениках "невинных молодых людей" и не станут им лгать.

Его ученики вернулись с ответом: "Для вас лучший вариант соединить несколько рабочих станций вместе и создать кластер; не надо ходить со своими большими запросами, к продавцам суперкомпьютеров". Шультен получил ответ! Ему придется собрать еще один параллельный компьютер, но это подразумевает значительно меньше работы, чем самодельный суперкомпьютер, который он построил в 1980-х годах.

Хорошо, что Шультен не был обольщен Connection Machine. Как раз в то время, когда Шультен принимал решение о том, как потратить свои деньги на оборудование, летом 1991 года The Wall Street Journal опубликовала статью о том, что некоторые суперкомпьютерные компании несправедливо субсидируются DARPA (Defense Advanced Research Projects Agency). DARPA было основано в Соединенных Штатах в 1958 году как ответ на неожиданный запуск спутника Советским Союзом. Миссия DARPA состояла в том, чтобы расширить технологии страны за пределы непосредственных потребностей военных. По сути, она действовала, чтобы убедиться, что неожиданность типа "Спутника" больше никогда не повторится DARPA было чистым продуктом холодной войны.

Оказывается, DARPA играла в фавориты с суперкомпьютерными компаниями, и компания, которая выпускала Connection Machine, была крупным получателем, на сумму 55 миллионов долларов. Когда эта история разразилась, это стало большим позором для администрации Буша (старшего). Кроме того, оказывается, что компания, создавшая Connection Machine, изобиловала бесхозяйственностью, паранойей и некомпетентностью. (Полный провал Connection Machine подробно описан Гэри Таубсом в этой статье). Без подпитки от DARPA компания начала терять деньги, и стало ясно, что их детище не выдержит производительности других конкурирующих суперкомпьютеров, популярных в то время. Клаус Шультен вовремя миновал минное поле холодной войны.

Но теперь Шультен оказался в немилости у суперкомпьютерного центра университета за то, что отклонил их предложение о сделке с Connection Machine. Но он ни разу не оглянулся. У него было видение, что компьютер может быть очень полезным инструментом в биологии, и это определило его решение построить самодельный кластер в 1993 году. "Это действительно важный инструмент для биологии. Я хотел сделать его полезным, и поэтому я должен был быть уверенным, что действительно служу своей области, а не своему эго." Это в основном обобщает принцип, которым он руководствовался при принятии решений. Ему нужен был лучший компьютер для вычислительной биофизики.

Изображение первого самодельного суперкомпьютера, построенного Клаусом Шультеном. Его группа будет продолжать строить еще много параллельных кластеров.Изображение первого самодельного суперкомпьютера, построенного Клаусом Шультеном. Его группа будет продолжать строить еще много параллельных кластеров.

Но компьютерный кластер, в который Шультен решил вложить деньги, был лишь частью уравнения. Его группе еще нужно было запустить молекулярную динамику на этом параллельном кластере. В 1994 году аспиранты Шультена работали над переносом существующих кодов молекулярной динамики, которые были разработаны предыдущими членами группы в конце 1980-х и начале 1990-х годов. Но студенты были так расстроены, глядя на эту массу неразборчивого кода. Они предложили написать новый код с нуля. Результатом стал NAMD, молекулярно-динамический софт, параллельный и написанный на современных C и C++. (Но можно долго рассказывать о презрении, которое они встретили за использование языков C и C++, а не почтенного ФОРТРАНА, которым пользовался в то время каждый уважающий себя ученый-вычислитель.) Имея собственные параллельные машины и новое программное обеспечение, специально разработанное для параллельных вычислений, группа была готова к началу нового тысячелетия.

Аквапорины, давно искомые водные каналы

В 1985 году Клаус Шультен стал свидетелем определения первой структуры мембранного белка. В то время он поставил перед собой цель имитировать белок в мембране. Примерно пятнадцать лет спустя, на рубеже веков, его группа наконец-то была готова ответить на вызов внедрить белок в мембрану. Кандидатом был скромный мембранный белок. Его функция состояла в том, чтобы просто впускать или выпускать воду из живой клетки. Но простота функций этого белка противоречит сложной сюжетной линии, которая раскрылась перед исследователями.

Вода необходима для всех живых клеток. Например, человек на 70% состоит из воды. Но до 1992 года ученые не могли объяснить, как некоторые клеточные мембраны позволяли проникать значительному количеству воды внутрь клетки. Ученые знали, что возможна некоторая диффузия через липидный бислой, но некоторые считали, что в таких клетках, как слюнные железы, почечные канальцы и эритроциты, проникновение происходит в большем масштабе, чем можно объяснить скромной диффузией. Так начались поиски водного канала.

Аквапорин встроенный в мембрану прокачивает воду внутрьАквапорин встроенный в мембрану прокачивает воду внутрь

Питер Агре называет свое открытие такого водного канала "чистой слепой удачей". В середине 1980-х годов Агре изучал антигены системы резус-фактора и наткнулся на небольшой мембранный белок, который он рассматривал просто как загрязнитель в своих препаратах. Лишь несколько лет спустя он начал задумываться о том, какова функция этого нового мембранного белка. Проведя множество различных исследований, чтобы узнать все, что можно о белке, он решил проконсультироваться с другими учеными, биохимиками и физиологами, чтобы получить подсказки. Один из его друзей, с которым он консультировался в Северной Каролине по дороге домой с летних каникул во Флориде, сказал, что это похоже на давно искомый водный канал.

Итак, Агре задумал эксперимент экспрессировать таинственный белок в систему, которая известна своей низкой водопроницаемостью. Он выбрал лягушачьи икринки, которые довольно стабильно живут в пресной воде. Когда он положил их в дистиллированную воду, заблаговременно одарив их таинственным белком, они полопались, как попкорн! Это означало, что таинственный белок позволял воде беспрепятственно проникать в систему, которая обычно счастливо живет в воде. Он нашел давно искомый водный канал, названный впоследствии "аквапорин". Свои результаты он опубликовал в апреле 1992 года.

К 2000 году были опубликованы две структуры аквапоринов, которые привлекли внимание Шультена: аквапорины человека и кишечной палочки. Для Шультена это была идеальная система, поскольку именно этот мембранный белок выполнял простую функцию: проводил воду. С точки зрения молекулярной динамики это позволяло Шультену роскошь наблюдать за мембранным белком в процессе его функционирования, что часто не было возможно для многих экспериментальных методов, используемых для характеристики белка. Затем в течение следующих семи лет в группе Шультена большое внимание уделялось моделированию аквапорина, что привело к появлению более пятнадцати публикаций на эту тему.

Пращи и стрелы

Ох уж эти теоретики! Они играют в свои игры и потом соглашаются с экспериментом. Они не говорят вам ничего нового. НЕТ, не в этом случае, оживляется Клаус Шультен, когда речь заходит о молекулярной динамике применительно к аквапоринам. Но он повторяет общую дилемму для ученых, занимающихся вычислительной биологией. Если ваши теоретические расчеты согласуются с экспериментом, то они не новы; а если ваши теоретические результаты не имеют экспериментального подтверждения, то они должно быть недействительны. Это практически безвыигрышная ситуация. Но работа Шультена с аквапорином действительно убедила его в том, что компьютер может вести себя почти как вычислительный микроскоп. Он часто может запечатлеть динамическую картину, где эксперимент слеп. В работе по аквапорину, которую проделал Шультен, он очень гордится одним открытием, в частности, тем, как эти каналы исключают протоны из потока воды. Но это была победа на пути к которой прошло множество битв.

К началу 2001 года Шультен уже подготовил к публикации ряд работ, посвященных механизму проводимости через аквапорины. Ранее было обнаружено, что некоторые аквапориновые каналы могут проводить глицерин так же, как и воду, и некоторые работы Шультена также были сосредоточены на этих так называемых акваглицеропоринах. Но Шультен понял, что ему нужны сотрудники-экспериментаторы, чтобы подкрепить открытия в области молекулярной динамики. Профессор замечает, что у него часто больше шансов быть опубликованным, если он объединился с экспериментальной группой одна из причуд жизни в вычислительной биологии. В это время у него был приглашенный аспирант из Дании Мортен Йенсен, работавший над акваглицеропорином кишечной палочки. Йенсен вместе с другим постдоком Эмадом Таджхоршидом проводили исследования и моделирование структуры аквапорина, определенной в 2000 году группой Роберта Страуда в Сан-Франциско. В конце концов Шультен убедил Страуда, что обе группы должны объединиться.

К тому времени Йенсен и Таджхоршид проводили моделирование на новой машине в Питтсбургском суперкомпьютерном центре под названием Terascale Computing System. В этот момент они действительно нуждались в ней, поскольку изучаемая ими система белок, мембрана, вода состояла из более чем 100 000 атомов. Но благодаря NAMD вычисления прошли плодотворно. И это в то время как многие исследователи могли себе позволить систему не более чем в 10 000 атомов.

Молекулы воды меняют свою ориентацию при прохождении по каналу аквапоринаМолекулы воды меняют свою ориентацию при прохождении по каналу аквапорина

Моделирование, проведенное Йенсеном и Таджхоршидом, подтвердило гипотезу (высказанную группой Есинори Фудзиеси в 2000 году) о том, что ориентация молекул воды может иметь какое-то отношение к магической способности аквапорина предотвращать протонную проводимость. В сущности, если есть цепочка молекул воды, протоны могут легко водородно связываться с водой, а затем прыгать от молекулы к молекуле, подобно тому, как продавщицы с чебуреками и мороженным переходят от вагона к вагону в движущемся поезде. По сути, моделирование показало, что вода движется по каналу кислородом вперед, а затем на полпути молекула воды переворачивается, и атомы кислорода смотрят против направления движения, как показано на рисунке. Боб Страуд предложил провести симуляции с отключенными ключевыми компонентами, которые, как подозревалось, могут быть ответственны за "сальто". Действительно, эти симуляции подтвердили наличие остаточных зарядов, ответственных за переворот.

Тетрамерная природа аквапорина.Тетрамерная природа аквапорина.

Но какой именно механизм отвечает за переворот и как это удерживает протоны от движения вместе с молекулами воды по каналу? По этому поводу велись длительные дискуссии. В конце концов ученые выяснили, что переворачиванию воды и удержанию протонов способствует симметрия белкового канала. В подразделе физики, электростатике, объект с противоположными зарядами на концах называется диполем. Два диполя вместе образуют квадруполь. Внутренняя часть аквапорина имела именно такое квадрупольное поле из-за своей симметрии. Наличие двух противоположных диполей эффективно мешает протонам прыгать вниз по цепочке воды, и поэтому протоны никогда не проходят через нее. Поскольку внутри клетки поддерживается почти постоянное напряжение, потеря протонов разрушит эту "клеточную батарею" и сделает клетку дисфункциональной. Кстати, позже Шультен обнаружил, что квадрупольное поле также участвует в проводимости глицерина. Глицерин длинная линейная молекула, и перевернуть молекулу, пока она находится в канале, было бы примерно так же легко, как развернуть на 180 градусов большой грузовой корабль, на самых узких участках Панамского канала (или Суэцого ха-ха). Молекула глицерина имеет три O-H группы, каждая из которых имеет свой собственный дипольный момент. И она, проходя через аквапорин, на самом деле просто вращает своими ответвлениями (которые содержат О-Н-группы) сверху вниз в квадрупольном поле. Таким образом глицерин, подстраивая свой дипольный момент к внутреннему квадруполярному полю движется по каналу.

С такими проработанными деталями на руках, полученными объединенными усилиями экспериментальной и теоретической команд, Шультен был поражен, когда редактор Science позвонил ему, чтобы объяснить, почему журнал отклоняет его статью. Судя по всему, рецензент раскатал ее в труху. И, чтобы подсыпать соль на рану, редактор сообщил, что вот-вот выйдет еще одна статья, тоже посвященная молекулярно-динамическому исследованию аквапорина. Команда Шультена была захвачена врасплох!

К счастью, Шультен смог договориться с редактором, и Science дала ему шанс на публикацию, но только в том случае, если команды проведут еще один тест, добавив указатели для определения положения воды. Итак, экспериментальная группа добавила меченные молекулы, а затем теоретики провели расчеты и получили очень хорошее согласие. Но это задержало публикацию на четыре месяца. В конце концов команды опубликовали свои результаты в апреле 2002 года.

Но на этом трудности Шультена не закончились. Более года спустя старый друг Шультена опубликовал статью, в которой утверждал, что теоретические результаты об исключении протонов, которые объяснил Шультен, были абсолютно неверны. "Работы теоретиков у людей всегда вызывают сомнения, рассказывает Шультен об одной из опасностей быть биологом-теоретиком. И как только другой теоретик утверждает, что работа ошибочна, то уже легко с этим согласиться. Наш прекрасный пример так и запоролся". Профессор решил ничего не делать и надеялся, что время покажет. Выжидательный подход в конце концов окупился сегодня его научная статья по аквапорину является довольно высоко цитируемой публикацией.

Для Мортена Йенсена, приглашенного датского студента в лаборатории Клауса Шультена, который провел так много исследований аквапоринов, работа, которую он сделал, имела глубокие последствия. Он рассказывает, что в лаборатории Шультена они обсуждали технологическую целесообразность использования аквапоринов для фильтрации воды. Когда Йенсен вернулся в Данию, он двинулся дальше, став соучредителем компании по разработке аквапоринов! Он набрал много ученых, сформировав консультативный совет, и не остановился на этом. "Мы уже получили стартовой финансирование, которое нужно только для того, чтобы сделать самые первые детские шаги, чтобы получить доказательство концепции" вспоминает Йенсен. С тех пор компания "Аквапорин" набирает обороты и сейчас работает над тем, чтобы расширить этот процесс для промышленных целей, таких как очистка сточных вод и, возможно, даже опреснение. В то время как опреснение в настоящее время требует больших затрат энергии для производства питьевой воды, датская компания разработала мембрану, полную белков, которая не требует высокого давления, чтобы протолкнуть воду через мембрану вместо этого ключом является осмос. Поскольку доступ к чистой питьевой воде является глобальной проблемой и по меньшей мере семь штатов в США страдают от сильной засухи, подобная технология имеет широкие последствия.

Для Клауса Шультена работа с аквапоринами стала настоящим открытием. Он всегда подозревал, что компьютер будет хорошим инструментом для мембранных белковых процессов, и в конце концов получил свое доказательство. Он увидел три ключевые области, где компьютер был полезен. Во-первых, теперь можно было достичь больших масштабов. Во-вторых, достаточно хорошо описана неупорядоченная мембранная среда. И, наконец, динамические процессы, происходящие в мембранных белках, такие как транспорт, могут быть воспроизведены, что особенно полезно, когда эксперимент не может дать динамических данных. Несмотря на многочисленные неудачи, у Шультена теперь был зеленый свет, чтобы двигаться вперед и описывать очередные мембранные белки.

Вход в новый класс белков: механочувствительные каналы

На заре тысячелетия были открыты новые структуры мембранных белков, что привлекло внимание Шультена. В то время как его группа активно участвовала в исследованиях аквапорина, он решил заняться другим классом мембранных белков: механочувствительными каналами. Этот объект исследования до сих пор ошеломляет Шультена своей загадочностью. Однако эта же загадка дает Шультену шанс на открытие, и он бросает на нее все, что имеет. Его любовь к открытиям, по общему признанию, является главным импульсом, который движет им как ученым.

Клетки прокариот и эукариот выработали способы реагирования на механический стресс. Например, бактериальная клетка любит, чтобы окружающая ее среда точно имитировала среду внутри нее. То есть концентрация ионов должна быть примерно одинаковой внутри и снаружи. Но что, если будет ливень, и вдруг снаружи бактериальной клетки окажется в основном вода, а внутри вода и растворенные вещества? Вода устремится в клетку, этот процесс называется осмосом. Но тогда клетка разбухнет и через некоторое время рискует взорваться из-за большого внутреннего давления. Это называется осмотическим шоком. Для избежания такого исхода и существуют механочувствительные каналы. Их три типа, и каждый из них срабатывает как своего рода предохранительный клапан, когда уровень тревоги достигает определенного значения, подобно тому, как аэропорты используют цветовой код для различных уровней тревоги в годы после 11 сентября. (К слову, в кишечной палочке был обнаружен четвертый канал, называемый MscK, но здесь будут обсуждаться только канонические три канала.)

MscL механочувствительный канал высокой проводимостиMscL механочувствительный канал высокой проводимости

Хотя эти механочувствительные каналы важны для клеток, подвергающихся осмотическому стрессу, они также вовлечены в осязание и слух. Кроме того, считается, что растительные клетки используют эти каналы для чувства гравитации и различения разницы между верхом и низом. В 1998 году экспериментаторы определили кристаллическую структуру механочувствительного канала, а к 2000 году Клаус Шультен уже проводил на нем моделирование. Этот канал называется MscL механочувствительный канал большой проводимости. Целью моделирования в группе Шультена было выяснить механизм стробирования или то, что заставляет канал открываться и закрываться. Эти симуляции, в основном выполненные аспирантом Джастином Гуллингсрудом, были первым моделированием молекулярной динамики в присутствии приложенного поверхностного натяжения. Это было введением Шультена в класс каналов, характеристики которых созрели для раскопок с помощью компьютерного моделирования.

Очарование воздушного шара

Как только Клаус Шультен увидел кристаллическую структуру MscS, он был чрезвычайно заинтригован он не мог поверить, насколько она прекрасна. Механочувствительный канал малой проводимости, или MscS, расположен в клеточной мембране бактерии образуя потенциальный проток. При нормальных условиях канал фактически закрыт, но при осмотическом напряжении, растягивающем клеточную мембрану, канал становится открытым. Открытие начинается на среднем уровне готовности; это второй из трех типов белков, которые открываются. По иронии судьбы, MscS это гораздо более крупный белок, чем MscL, канал большой проводимости. На самом деле большой размер MscS это одна из тех вещей, которые одновременно завораживают и сбивают с толку Шультена. Большая часть этой конструкции (около 65%) находится вне клеточной мембраны, в цитоплазме. Как показано на рисунке, MscS во многом напоминает фонарь, нижний сегмент которого находится вне липидного бислоя и содержит семь отверстий-окон. Эта нижняя часть, или цитоплазматический домен, была прозвана в группе Шультена "воздушным шаром". Когда Шультен начал изучать MscS в 2003 году, он понял, поскольку все обращали внимание на часть белка, которая находилась в липидах (трансмембранный домен), что природа и функция баллона остались загадкой. И он хотел ее разгадать. Многие считают, что воздушный шар это своего рода фильтр. К сожалению, первые набеги Шультена на MscS, казалось, давали больше вопросов, чем ответов.

Поскольку Шультен уже изучал механочувствительные каналы (в форме MscL), когда новый аспирант Маркос Сотомайр присоединился к группе примерно в 2003 году, и впервые появилась кристаллическая структура MscS, казалось естественным для Шультена продолжить эту тему. Сотомайр говорит, что это было для него хорошей проблемой для введения в биофизику белков.

MscS с "баллоном". В шаре видны два окна.MscS с "баллоном". В шаре видны два окна.

К счастью, Сотомайр уже был хорошо знаком с молекулярной динамикой: он написал свою собственную программу для работы над темой, связанной с физикой фазового перехода в газах. Он проводил эту работу, когда был студентом-физиком в Чили. Сотомайр рассказывает, что одной из причин, по которой он выбрал Университет Иллинойса в Урбана-Шампейне для аспирантуры, был широкий спектр областей, охватываемых более чем шестьюдесятью профессорами физического факультета. Во время поиска научрука, Сотомайр увидел лекцию Клауса Шультена о том, как птицы используют квантовую механику для ориентации в магнитном поле земли. Это настолько заинтриговало Сотомайра, что он обратился к профессору. В то время как студент был немного обеспокоен тем, что у него не было опыта в биофизике, Шультен был в восторге от того, что у него была подготовка по молекулярной динамике "мне нужны люди, которые знают, как делать симуляции, а биологию вы узнаете по пути".

Когда Сотомайр приступил к моделированию MscS, у всех, включая его самого, сложилось впечатление, что кристаллическая структура 2002 года была зафиксирована в открытом состоянии. В принципе, когда вода + ионы проходят через пору, это представляет собой открытое состояние; когда ионы не могут пройти, это представляет собой закрытое состояние. И существует третий вариант неактивное состояние. Раннее моделирование, проведенное Сотомайром и Шультеном на MscS, показало асимметричное закрытие поры. Здесь наши герои получили свою первую подсказку о том, что кристаллическая структура 2002 года не может быть одним из открытых состояний. Наряду с обнаружением асимметричного замыкания они также обнаружили, что липиды взаимодействуют с белком весьма интенсивно, поэтому отношения между мембраной и белком были сложными. На самом деле, по рассказам Сотомайра, сложное липидно-мембранное взаимодействие "фактически определило это очень динамичное движение трансмембранного домена, которое привело к асимметричному замыканию." В сущности, Сотомайр и Шультен видели большую деформацию липидной мембраны вокруг белка. Также было обнаружено, что поры в трансмембранном домене несколько гидрофобны, то есть вода будет проникать в поры, но только периодически.

На следующем этапе работы над MscS Сотомайр и Шультен решили воспроизвести некоторые экспериментальные работы. Поскольку MscS это канал, по которому протекают ионы, экспериментаторы измерили проводимость, и эти измерения можно было воспроизвести в симуляциях. Они объединились с другим ученым университета, Умберто Равайоли и его постдоком Труди ван дер Страатеном. Ученые обнаружили, что кристаллическая структура, которую они использовали, действительно не была открытой. Таким образом, теоретики попытались воспроизвести возможные открытые формы канала в своих расчетах. "И я думаю, что это было действительно важно, потому что это направляло будущие эксперименты на попытки найти более широкий открытый канал", объясняет Сотомайр.

Сотомайр согласен с тем, что эти первые два исследования, в результате которых были подготовлены две статьи, казалось, вызвали больше вопросов, чем дали ответов. "Это было правильное исследование, говорит Сотомайр о своих первых двух публикациях, и мы просто смотрели на что-то новое, а затем находили вещи, о которых мы точно не знали, что они означают." Во-первых, эти первые две работы предполагали, что кристаллическая структура 2002 года была либо неактивным состоянием, либо закрытым состоянием, но теоретики не были уверены, какое именно. В то же время было неясно, какова структура или конформация открытого канала. В группе Шультена, в работе, выполненной тогдашним постдоком Алеком Аксиментьевым, был разработан метод изучения электростатического потенциала через пору. Сотомайр рассказывает, что когда он рассчитывал электростатический потенциал для MscS, он вышел странным появилась еще одна загадка. Также без ответа остался вопрос о том, как именно функционирует механизм стробирования, каковы все переменные, способствующие переходу из открытого состояния в закрытое. А потом был воздушный шар. Какова его функция и как работают семь боковых окон? Играл ли воздушный шар большую роль в стробировании?

Есть много причин, почему белковые MscS так сложны как для экспериментаторов, так и для ученых-вычислителей. Это большой белок, в общей сложности более тысячи аминокислот. Он состоит из семи субъединиц. Как упоминает Сотомайр, с таким большим размером трудно точно определить, какие участки наиболее важны для функции белка. И он действует в клетке в миллисекундном масштабе времени, а современное моделирование позволяет охватить сотни наносекунд. Еще в 2003-2005 годах, когда Сотомайр проводил моделирование, он достиг предела размера в 224 000 атомов и мог моделировать только несколько наносекунд,

Различные состояния мембранных белков MscSРазличные состояния мембранных белков MscS

Когда Сотомайр и Шультен опубликовали свои первые две работы, у них было большое подозрение, что кристаллическая структура, которую они использовали, вероятно, не была открытым состоянием, но они нуждались в более глубоком анализе, чтобы исследовать эту дилемму. На самом деле это яркий пример в карьере Шультена, когда компьютер может помочь в определении различных состояний мембранных белковых каналов, поскольку некоторые из этих состояний не обязательно поддаются кристаллизации. Но чтобы получить наилучшие результаты, придется объединиться с экспериментальной группой. Так началось плодотворное сотрудничество с другом Эдуардо Перозо, специалистом в области молекулярной физиологии и биологической физики, который сосредоточил свои исследования на динамике ионных каналов. Студентка Перозо, Валерия Васкес, в то время работала над MscS, и поэтому она и Сотомайр часто кооперировались. Сотомайр должен был убедиться, что моделирование, которое он проводил, совпадает с экспериментальными условиями, которые Васкес создавала в лаборатории. Затем обе команды приступили к трудной работе по расшифровке точного значения кристаллической структуры MscS. Фактически, в своей первой совместной публикации они написали:

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

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

Как оказалось, сотрудничество Шультена и Перозо оказалось плодотворным. Например, в то время как Васкес измеряла проводимость через канал MscS, Сотомайр был в состоянии воспроизвести измерения в моделировании. Объединенные команды также смогли объяснить стробирующее поведение MscS, что, по сути, привело к публикации в Science. Но самое главное, сотрудничество подтвердило, что первоначальная кристаллическая структура была не открытой, а скорее неактивной или закрытой. Сотомайр и Шультен смогли учесть экспериментальные данные из лаборатории Перозо и создать некоторые возможные модели открытого состояния для MscS, что является классическим использованием компьютера при разработке состояний мембранных белков. На самом деле, Сотомайр считает, что экспериментальные данные действительно были "ключом к поиску открытых состояний."

Как оказалось, как раз в то время, когда сотрудники опубликовали свою научную работу 2008 года о стробирующем механизме MscS, в Science также появилась статья, представляющая чрезвычайный интерес для команд Шультена-Перозо. Экспериментальная группа, возглавляемая Яном Бутом и Джеймсом Нейсмитом, фактически определила кристаллическую структуру MscS с порой, намного более широкой, чем исходная кристаллическая структура 2002 года. Сотомайр говорит, что теоретические предсказания открытого состояния, предложенные им и Шультеном, которые в значительной степени опирались на эксперимент, очень напоминают открытое состояние структуры 2008 года. Но, как обычно с MscS , все детали открытого состояния все еще остаются предметом споров и некоторой путаницы.

Методологическая разработка для Энигмы

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

Прямо сейчас Шультен и его аспирант Иван Тео находятся на ранней стадии проекта по определению функции воздушного шара, и они очень взволнованы этим. "То, что он делает," говорит Шультен о работе Тео, "очень интересно не только с точки зрения MscS, но и с методологической точки зрения." Тео и Шультен объединили усилия, чтобы раскрыть секреты воздушного шара, разрабатывая общий метод для применения к MscS, но который достаточно общий, чтобы его можно было применить к другим мембранным процессам. И их работа является ярким примером того, как компьютер может захватить мириады деталей биологии.

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

Продолжение следует...

Подробнее..

Перевод Суперкомпьютеры и клеточные мембраны (заключительная часть)

07.05.2021 10:12:10 | Автор: admin

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

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


Предыдущие части: первая, вторая, третья.

Оглавление

  1. Транслоконы: дамбы Теночтитлана

  2. YidC: продолжение для белково-проводящего канала

  3. BAR-домены и мембранная скульптура

  4. Понимание нервной системы

  5. Грандиозный финал: серия взаимосвязанных процессов

Транслоконы: дамбы Теночтитлана

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

Теночтитлан, с фрески Диего Риверы в Национальном Дворце Мексики. В середине изображения находится дамба, ведущая из города-государстваТеночтитлан, с фрески Диего Риверы в Национальном Дворце Мексики. В середине изображения находится дамба, ведущая из города-государства

У жителей также были формы валюты (например, какао-бобы), которые гомологичны АТФ. Подобно тому, как аквапорины контролируют поток воды в клетке, ацтекский город имел каменный акведук для питьевой воды. Но как основные блага попадают в клетку и выходят из нее? В Теночтитлане было три дамбы. В клетке самым фундаментальным благом, вероятно, являются ее белки. А проходы в клетке, которые позволяют проходить этим критическим молекулам, белкам, называются транслоконами. Транслокон сам по себе является белком, который живет в мембране, и выполняет функцию проводящего канала. Когда Клаус Шультен начал изучать транслоконы в 2005 году, он не мог представить, что его ждет. Это как не знать о всех сокровищах затерянного города, пока они не откопаны, и Шультену придется раскопать некоторые древние исследования, чтобы изучить новую систему, которую он вскоре стал считать сокровищем.

Когда Джеймс Джей Си Гумбарт в 2003 году поступил в аспирантуру в Урбана-Шампейн, он ничего не знал о биофизике и предполагал, что в конечном итоге изучит физику конденсированных сред, что было сильной стороной университета. "Но потом, на первом семестре, я пошел на биофизику и был очень заинтригован ею", вспоминает Гумбарт. Это было лишь вопросом времени, когда студент начнет работать с Шультеном, поскольку вычислительная направленность также импонировала ему.

Летом 2004 года, когда Гумбарт присоединился к группе, Шультен показал ему статью о канале белковой проводимости, которая, по мнению Шультена, была интересной и, возможно, могла бы стать плодотворным проектом. Вероятно, к восторгу Шультена, Гумбарт действительно обратил на это внимание, проявил большую инициативу и продолжал производить впечатляющий объем работы.

Джей Си Гумбарт, Клаус Шультен и рибосома, прикрепленная к проводящему белок каналу.Джей Си Гумбарт, Клаус Шультен и рибосома, прикрепленная к проводящему белок каналу.

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

Примерно в то время, когда Шультен изучал транслоконы, его научные интересы внезапно расширились и включили в себя увлекательную и массивную молекулярную машину рибосому, биологическую фабрику, которая собирает белки клетки. В 2004 году один экспериментатор обратился к Шультену на конференции и попросил его помочь решить дилемму, связанную со структурами рибосом, полученными двумя различными типами микроскопов. Один тип микроскопа включает в себя рентгеновскую кристаллографию, которая выдает структуру высокого разрешения; другой тип электронная микроскопия, которая генерирует структуры низкого разрешения. Этот экспериментатор (Иоахим Франк) получал структуры рибосом с помощью криоэлектронной микроскопии и нуждался в помощи, включающей данные высокого разрешения (от рентгеновской кристаллографии, сгенерированной другими) в его карты. По аналогии рассмотрим две фотографии: на одной птица сидит на проводе, а на другой птица в полете. Изображение птицы в полете захватывает птицу в процессе полета (ее функция), но имеет низкое разрешение. Однако представьте себе, что первый снимок может запечатлеть мельчайшие детали, вплоть до каждого пера. Если бы каким-то образом можно было использовать все мелкие детали неподвижной птицы и применить их к птице в полете, можно было бы объяснить функцию птицы полет. Поэтому задача Шультена состояла в том, чтобы взять данные рентгеновской кристаллографии и поместить их в карты криоэлектронной микроскопии; эти карты электронной микроскопии были действительно способны захватить рибосому в процессе выполнения ее функции, в отличие от рентгеновской микроскопии.

Шультен оценил свои силы и согласился на проект. Его группа работала над большими системами в течение десятилетий, у него были собственные программные продукты в качестве инструментов, и у него была группа студентов, одаренных в вычислительном искусстве. Результатом стал метод, названный MDFF или молекулярная динамика с адаптивной подгонкой, и он дебютировал в 2008 году. Более подробную информацию о генезисе MDFF можно найти в этой истории. Достаточно сказать, что она открыла Шультену захватывающий мир рибосомы.

Но как рибосома связана с проводящим белок каналом? Как уже отмечалось, рибосомы производят белки, которые в процессе построения выходят из своего туннеля, а затем должны свернуться и отправиться в свой дрейфующий путь. Один из вариантов заключается в том, что зарождающийся белок может жить в цитозоле клетки. Или его судьбой может быть путешествие через мембрану наружу. Третий вариант зарождающийся белок может быть предназначен для жизни внутри мембраны. Для этих двух последних вариантов белок-проводящий канал направляет зарождающийся белок к цели. В 2007 году Шультену и Гумбарту стало известно о группе, которая опубликовала карту рибосомы в комплексе с ее белково-проводящим каналом полученную с помощью криоэлектронной микроскопии. (Этот канал известен как SecY). Это интригующая система сразу две макромолекулы, каждая из которых влияла на другую. Более подробно о функции рибосомы можно было бы узнать, если бы ее изучали в тандеме с проводящим белок каналом, который на нее влиял. И Шультен и Гумбарт поняли, что они могут генерировать структуры высокого разрешения каждой части дуэта, чтобы в конечном итоге скорректировать карту электронной микроскопии. Это казалось идеальной работой для MDFF. Гумбарт даже отправился в Бостон, чтобы получить карту электронной микроскопии от экспериментаторов, создавших комплекс.

Это было одно из первых применений нового метода MDFF для группы Шультена. Гумбарт рассказывает, что группа все еще отлаживала метод, когда он работал над этим комплексом рибосомы плюс канал в январе 2008 года, но это было очень захватывающее время, особенно приятно было видеть, как структуры почти волшебным образом вписываются в электронную плотность. MDFF скомпоновал комплекс из 2,7 миллиона атомов, который включал рибосому, белок-проводящий канал, мембрану и воду самый большой на тот момент в группе Шультена. Исследователи действительно могли видеть, как рибосома вызывала небольшую дестабилизацию в пробке канала. Это был очень успешный проект, и он был лишь прелюдией к тому, что должно было произойти.

В процессе работы над MDFF Шультен скооперировался с экспериментатором в Мюнхене, биохимиком по имени Роланд Бекман, который сумел выловить активирующие системы с помощью своего электронного микроскопа. Бекман уже работал над электронной микроскопией рибосомы, соединенной с проводящим белок каналом, когда он и Шультен объединили усилия. На самом деле, проект Шультена и Бекмана принес обеим командам научную публикацию, в которой MDFF был использован Гумбартом для карты канала связанного с рибосомой.

Клаус Шультен навсегда запомнит свой первый визит к Бекману в Мюнхене. Это было в июле 2008 года, и во время встречи Бекман вывел Шультена на улицу, чтобы показать ему плакат на стене с изображением электронной микроскопии. И Шультен чуть не упал в обморок от восторга, увидев кое-что знакомое. Он сразу же узнал на карте нанодиск, объект, смоделированный в университете Иллинойса. Шультен помогал открывателю нанодиска визуализировать его, так что это было то, над чем Шультен работал в начале своей карьеры. В принципе, открыватель нанодиска, Стив Слигар, хотел изготовить наноразмерный кусок мембраны, чтобы удерживать мембранные белки, поскольку, как уже отмечалось, работа с мембранными белками вне мембраны представляла экспериментаторам много проблем. Итак, Слигар дал миру рецепт создания пучка липидов, удерживаемых вместе двумя белками-каркасами маленький кусочек мембраны, на который можно посадить белок!

Рибосома и нанодиск - сложная, но увлекательная системаРибосома и нанодиск - сложная, но увлекательная система

Новшество Роланда Бекмана состояло в том, что он взял рибосому, в которой зарождающийся белок еще не полностью вышел из туннеля рибосомы, и ввел в систему нанодиск затем зарождающийся белок змеился в канале в нанодиске. В сущности, Бекман поймал моментальный снимок рибосомы с зарождающимся белком, входящим в канал в окружении мембраны. И Клаус Шультен сразу узнал нанодиск. Это о чем-то говорит, потому что, когда Шультен помогал Слигару визуализировать нанодиск еще в 2005 году, он получил лишь косвенную информацию о его форме. В сознании Шультена нанодиск выглядел очень похожим на то, что он видел на стене у Бекмана в тот летний день 2008 года, так что это было подтверждением того образа, который Шультен представил научному сообществу.

Бекман хотел, чтобы опыт Шультена в MDFF позволил создать структуру системы рибосома-нанодиск с высоким разрешением. Гумбарту казалось естественным работать над этим проектом, поскольку он включал в себя проводящий белок канал. Однако, по мнению Шультена, это был сложный проект. А ведь аспирант только начал учиться разбираться в картах электронной микроскопии. "Мне пришлось смотреть на этот канал часами, прежде чем я смог действительно начать мысленно видеть белок, выходящий из него", рассказывает Гумбарт. Но, к счастью, с ним работали еще два студента-пионера MDFF, чтобы помочь ему разобраться в картах.

Гумбарт и Шультен в конце концов смогли разобраться в сверхсложной, но новаторской системе, которую Роланд Бекман запечатлел с помощью электронного микроскопа. Для Гумбарта это было кульминацией всей его диссертационной работы. Одна из целей, которую они с Шультеном поставили перед собой, состояла в том, чтобы увидеть, как образуются мембранные белки. Для Шультена этот проект с рибосомой и нанодиском был вершиной его карьеры. Он считает, что этот проект подчеркивает, как вычислительная техника (через MDFF) может стать ключевым партнером для экспериментов. Но самое главное, он никогда не ожидал, что снова вернется к нанодиску, и он никогда не думал, что это будет что-то настолько инновационное.

YidC: продолжение для белково-проводящего канала

Скорее всего, проект провалится. Несмотря на это предупреждение Клауса Шультена, постдок Абхишек Сингхарой решил уделить ему все свое внимание. Летом 2013 года Шультен посетил Мюнхен и вернулся с несколькими новыми проектами, один из которых был совместным с группой Роланда Бекмана. Несмотря на мрачное напутствие, Шультен, вероятно, заинтересовался этой темой, потому что она была продолжением его многолетней работы над белковыми проводящими каналами. Для полного удовлетворения рибосомы, зарождающегося белка и транслокона было недостаточно, поскольку было известно, что дополнительный мембранный белок, называемый YidC, иногда помогает каноническому белковому проводящему каналу. Это было мучительно для Шультена, который всегда думал о системах с несколькими белками. Но он полностью отдавал себе отчет в шансах на успех, когда ставил перед Сингхарой задачу: экспериментаторам нужна помощь вычислителей, завершающих работу над структурой YidC.

Видите ли, для YidC не существовало никакой кристаллической структуры. И не было похоже, что она появится в ближайшие годы или десятилетия. Но Роланд Бекман считал, что структура может быть построена даже без дифракции рентгеновских лучей. Экспериментаторы в лаборатории Бекмана, после того как не смогли получить исходную структуру для YidC, обратились к вычислительной группе в своем родном университете (Мюнхенский университет Людвига Максимилиана), известной своими программными инструментами, предназначенными для предсказания структуры белка по последовательности. И на Сингхароя легла следующая задача работать рука об руку с группой Бекмана над усовершенствованием структуры. Сингхарой использовал бы все инструменты из арсенала группы Клауса Шультена, а экспериментаторы Бекмана, и в первую очередь аспирант Стефан Уиклз, сосредоточились бы на электронной микроскопии и биохимии. Это был прекрасный пример работы компьютера в тандеме с экспериментом. И это было путешествие, полное драм.

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

YidC, которая облегчает введение белка в мембрану.YidC, которая облегчает введение белка в мембрану.

Имея под рукой "входную" структуру, Сингхарой поместил белок YidC в мембрану и начал проводить моделирование молекулярной динамики. "Но это не было обычным моделированием молекулярной динамики", рассказывает Сингхарой. "Это была динамика ограниченная полученными ранее ковариационными измерениями". После того, как структура была запущена с ограничениями, затем без них, чтобы уравновесить и отобразить расслабленную структуру, пришло время отправить результаты сотруднику-экспериментатору, который проверяет, какие взаимодействия делают эту расслабленную структуру стабильной. Чтобы провести всю эту биохимию, Уиклз менял компоненты и тем самым устранял взаимодействия, которые предполагались как стабилизирующие. И потом весь процесс сызнова молекулярная динамика с ограничениями и без них, затем биохимические эксперименты. Сингхарой рассказывает, что таким образом они прошли через множество итераций. Все с целью найти наиболее устойчивую структуру.

После этих утомительных повторений Уиклз нашел способ удостовериться, что очередная структура Сингхароя была истинным энергетическим минимумом, а не каким-то локальным, в котором задержался YidC, потому что время моделирования было недостаточно длительным. Уиклз провел химический эксперимент по сшиванию. В принципе, если два остатка (парочка определенных аминокислот) расположены близко, то после добавления к ним группы серы, они могут образовывать дисульфидную связь, и это может быть сравнительно просто обнаружено, таким образом доказывая, что два рассматриваемых остатка находятся рядом. Чтобы проверить структуру Сингхароя, Уиклз ввел в YidC зарождающийся белок и перекрестно связал некоторые их остатки. И о чудо, они были связаны! Это был самый сильный тест, который показал, что многие итерации окупились, и структура была значимой.

И в разгар этой тяжелой работы, Уиклз узнал немыслимое: их нагоняла другая группа! Вот-вот должна была быть опубликована кристаллическая структура бактериального YidC. Уиклз был на научной конференции, где у него был плакат о его работе над YidC, когда он увидел кристаллическую структуру на другом плакате. Если и было что-то хорошее, так это то, что Уиклз заметил, что его смоделированная структура очень похожа на кристаллическую структуру, которую он видел на конференции.

Это вызвало волну лихорадочных приготовлений к подготовке публикации. Три команды, в принципе, уже были готовы представить доклад, когда Уиклз увидел кристаллическую структуру на конференции. Наконец, после долгих испытаний и невзгод, статья была принята в июле 2014 года и вышла только через три месяца после того, как кристаллическая структура была опубликована командой из Японии. Структура, основанная на модели, которую опубликовала триада команд, была не единственным значительным результатом работы. Благодаря обширной молекулярной динамике, которой руководил Сингхарой, он начал раскрывать некоторые детали того, как на самом деле работает YidC, поскольку он мог видеть его динамически в своих симуляциях. Во-первых, моделирование показало, что дно YidC образует уплотнение. Это имеет смысл, поскольку известно, что YidC является инсертазой, то есть она берет зарождающуюся цепь из рибосомы и направляет эту цепь в сторону ее последнего пристанища в мембране. В процессе, ни один канал не открывается внизу. Вместо этого новообразованный белок проходит через ворота в боку YidC. Моделирование показало, что одна из спиралей, образующих ворота, действует очень похоже на затворку или гибкий закрылок, позволяющий зарождающемуся белку змеиться в мембрану.

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

BAR-домены и мембранная скульптура

Новая аспирантка Ин Инь хорошо помнит тот день в 2006 году, когда Клаус Шультен позвал ее в свой кабинет и достал из кармана смятый клочок бумаги. Он протянул ей листок и сказал: "Мы должны поработать над этим", на бумаге она увидела только три буквы: B...А...R. Не так уж много, для уверенного старта. Но в течение следующих нескольких лет эти три литеры приведут к проекту, в котором будут задействованы все инструменты из арсенала Schulten group. И созвездие факторов сошлось вместе, чтобы произвести прекрасное применение молекулярной динамики. Этот проект, на самом деле, является непревзойденным примером того, как вычислительная биология может дать представление о динамических процессах, связанных с мембранными белками. Ниже мы увидим, как такая клеточная активность как ваяние мембран, была полностью раскрыта от начала до конца.

Что же такое мембранная скульптура? Эукариотическая клетка это сложное строение с множеством отсеков, отгороженных мембранами. Эти мембраны бывают самых разных форм плоские, изогнутые, волнистые или правильной формы. Есть также везикулы в клетке, например, в аппарате Гольджи (amazon.com клеточного мира) белки упаковываются внутри мембран (один тип везикул) и отправляются в свой конечный дом, либо в цитоплазму, либо, возможно, высвобождаются наружу клетки. Дело в том, что внутри клетки существует множество мембранных систем, и каждая из них имеет уникальную форму. Но что именно отвечает за появление столь многих различных форм? Часто для "лепки" мембраны в ее окончательный вид используются белки.

Когда Клаус Шультен дал Ин Инь этот листок бумаги с тремя магическими буквами, он на самом деле подразумевал семейство белков, называемых BAR-доменами. Известно, что эти белки изгибают мембраны (во время эндоцитоза и экзоцитоза). На самом деле существует три члена семейства (N-BAR, F-BAR, I-BAR), и в центре внимания Инь вскоре оказались белки N-BAR, для которых уже была известна структура. По мере того как она начинала исследовать белок, она становилась все более и более заинтересованной и все более и более убежденной, что это тема, достойная изучения. И Инь начала обсуждать N-BAR-домены со своим женихом Антоном Архиповым, тоже аспирантом в группе Шультена. Супруги поняли, что могут сотрудничать и тем самым промоделировать то, что раньше казалось невозможным. Даже Шультен не ожидал, что такое моделирование осуществимо. О чем эта парочка догадалась?

Инь и Архипов познакомились в первые же дни учебы в аспирантуре Иллинойского университета в августе 2004 года. Они пересеклись на квалификационном экзамене и поняли, что оба новички на физическом факультете, Инь из Китая и Архипов из России. Архипов, увлеченный нейробиологией, присоединился к группе, поскольку он узнал о прошлых работах Шультена в области нейробиологии, а затем увлекся идеей вычислительной биофизики, над которой Шультен в настоящее время работал. Инь начала работать в лаборатории физики высоких энергий, но поняла, что хочет заниматься чем-то связанным с жизнью, и Архипов предложил ей присоединиться к группе Шультена.

Сравнение полноатомного представления по сравнению с крупнозернистым на одном белке BAR-доменаСравнение полноатомного представления по сравнению с крупнозернистым на одном белке BAR-домена

У Архипова было много разнообразных проектов в группе Шультена, и над одним из них он работал вместе с аспирантом Питером Фреддолино, а именно над реализацией крупнозернистой молекулярной динамики в NAMD. В этом методе группы атомов объединялись в эдакие "бусинки". Например, одна бусина может представлять 10 атомов. Хотя это может означать некоторую потерю мельчайших атомных деталей, это, по существу, равносильно возможности управлять молекулярной динамикой в чрезвычайно больших системах и в течение чрезвычайно длительного времени. Архипов занимался двумя видами крупнозернистой молекулярной динамики, одна из которых была основана на остатках, а другая на форме. Первый имеет разрешение одного остатка (около 10 атомов на шарик), а второй использует группу шариков для представления целого белка (около 150 атомов на шарик).

Архипов фактически использовал крупнозернистую молекулярную динамику на вирусных капсидах в качестве приложения метода. Также в это время аспирантка Эми Ши использовала крупнозернистую молекулярную динамику для упомянутой выше системы нанодисков. "Поскольку Инь работала над BAR-доменами, я, конечно, обсуждал это с ней просто из любопытства", рассказывает Архипов. "И мы вместе решили, что это может быть хорошим применением крупнозернистой динамики." Супруги рассказывают, что Шультен очень поддержал их предложение о сотрудничестве. К слову, Шультен сотрудничает со своей женой, химиком, уже около сорока лет и понимает как выгоды, так и потенциальные профессиональные опасности сотрудничества мужа и жены в науке. Инь и Архипов поженились в 2008 году.

Кривизна мембраны создается несколькими BAR-доменамиКривизна мембраны создается несколькими BAR-доменами

В первой работе Инь и Архипов поместили один доменный белок, который имеет форму более или менее похожую на банан, на участок мембраны и наблюдали, как молекулярная динамика показала, что мембрана рядом с белком изгибается, чтобы соответствовать кривой серповидного белка. На следующем шаге, они поместили шесть BAR-белков в линию на участке мембраны. Единственная причина, по которой они могли просто попытаться совершить такой подвиг, была новая крупнозернистая функция NAMD. Без упрощенного представления этот финт потребовал бы обсчитать несколько миллионов атомов, что не позволило бы моделировать достаточно долго, чтобы увидеть, согнули ли шесть белков мембрану. Когда шесть стержневых белков образовали волнистую форму, они попытались сложить шесть белков в виде решетки. Вуаля! Они наткнулись на конфигурацию, которая произвела красивую кривизну. Это была решетка из BAR-белков.

Шультен рассказывает, что, когда они попытались опубликовать свои результаты с решеткой BAR-доменов, у них возникли проблемы. Эксперты отмечали, что для клеток не характерны высокие концентрации BAR-белков, поэтому моделирование казалось нереалистичным. "Но потом кому-то удалось сделать электронные микрофотографии мембран вместе с BAR-доменами", отмечает Шультен. "А потом они увидели, что те образуют именно ту решетку, о которой мы говорили!" Так что теперь казалось вероятным, что многие BAR-домены работали вместе, чтобы лепить мембраны. Вскоре после этого работа была принята.

В качестве грандиозного финала Шультен, Инь и Архипов решили действительно проверить новый метод крупнозернистую молекулярную динамику, и посмотреть, смогут ли они сформировать трубку из плоского участка мембраны. "Психологически, и просто как демонстрация силы метода, ничто не может быть лучше, чем фактически показать формирование трубки целиком", отмечает Архипов. Поэтому Инь поместила сеть из N-BAR-белков на кусок мембраны в 200 нанометров в квадрате, передала свою работу суперкомпьютеру и стала ждать. И ждала она почти 200 дней. Однако Инь рассказывает, что она проверяла прогресс несколько раз в день и часто видела одно и то же изо дня в день трубка не закручивалась. Шультен велел ей набраться терпения. И вот однажды образовалась идеально закрытая трубка! Все трое ученых были в восторге. Крупнозернистая молекулярная динамика могла дать масштабный обзор процесса лепки мембраны. Как подытоживает Архипов: "Я был очень рад видеть это, потому что изобрести метод легко, но сделать его полезным сложнее."

Шультен продолжает свою работу с BAR-доменами и по сей день. В настоящее время его аспирант Ханг Юй работает над F-BAR. Вдохновленные полученным в 2008 году криоэлектронным микроскопическим изображением решетки F-BAR-доменов на трубке, Юй и Шультен приступили к ответу на ряд вопросов, к примеру: как работает F-BAR-домен в клетке, лепит ли он мембрану подобно N-BAR-белкам? Почему клетка использует такую специфическую концентрацию F-BAR-доменов?

Юй решил посмотреть, как F-BAR-домен изгибает кусок мембраны. Когда он запустил свою симуляцию, ничего не произошло, мембрана просто осталась плоской. Он провел несколько симуляций и так ничего и не увидел. Но он использовал те же времена, что и его предшественники, Инь и Архипов. Юй, наконец, решил позволить моделированию работать очень долго, и, о чудо, он наконец увидел, как F-BAR-домен изгибает мембрану. Юй обнаружил, что F-BAR-домены менее жестки, чем N-BAR-домены, но в то же время производят меньшую кривизну.

Затем Юй попытался найти, какая оптимальная плотность белков F-BAR-домена дает наибольшую кривизну. На кусочек мембраны размером 1000 квадратных нанометров он положил 5, 8, 10, либо 16 димеров. Он увидел, что расположение 10 димеров создает наилучшую кривизну на плоской мембране. Теперь Юй был вооружен особой решеткой, которая производит наибольшую кривизну. Он был готов посмотреть, будут ли F-BAR-домены производить трубочку из плоского куска мембраны. Юй сообщает, что он и Шультен были в восторге, когда увидели, что трубка полностью сформирована.

Решетка из BAR-доменов сминает мембрану в трубкуРешетка из BAR-доменов сминает мембрану в трубку

Во всех работах по BAR-доменам Шультен утверждает, что то, что на самом деле создает кривизну мембраны, называется механизмом строительных лесов. В этом механизме есть положительно заряженные остатки на нижней стороне банановидного BAR-домена, которые притягивают отрицательные липидные головки. Это притяжение смещает липиды к BAR области и таким образом искривляет мембрану. Но многие все еще не определилось с механизмом лепки, поскольку некоторые исследователи предпочитают другой метод, называемый механизмом вставки. В этом механизме остатки BAR-домена фактически вставляются в мембрану и зажимают ее, заставляя изгибаться. Шультен признает, что иногда в клетке некоторые белки могут использовать механизм вставки для изгиба мембран. "Потому что если заглянуть в клетку, поясняет Шультен, то можно увидеть очень много разных форм. Я имею в виду так много, что вы никогда бы не подумали, что все они появляются в с применением только одного трюка."

Понимание нервной системы

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

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

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

Примерно в то же время, у него также появился новый ученик, Фатима Халили-Арагхи, которая взялась изучать эту тему. Первоначально заинтересовавшись физикой высоких энергий, при поступлении в аспирантуру она вскоре поняла, что физика в группе Шультена была такой же увлекательной, как и в теории струн. Она немедленно начала исследовать калиевые каналы. Ученые знали, что это канал, и он открывается и закрывается (эффект стробирования) в зависимости от изменения напряжения. "Когда я начала изучать литературу, я поняла, что мы не знаем, как они чувствуют напряжение. Мы не знали, что происходит внутри", вспоминает Халили-Арагхи. Ее задача состояла в том, чтобы узнать больше и прояснить механизм.

Калиевый канал встроенный в мембрану активируется при возникновении напряженияКалиевый канал встроенный в мембрану активируется при возникновении напряжения

Она и Шультен впервые попытались понять проникновение ионов калия через поры. В белке есть пять аминокислот, которые образуют самый узкий сегмент канала и которые не пропускают ничего, кроме калия, поэтому он известен как селективный фильтр. Многое было постулировано о том, как ионы проходят через селективный фильтр, и одна идея заключалась в том, что калий не течет непрерывно. Когда Халили-Арагхи начала имитировать калий идущий потоком, канал не проводил. Это было странно, потому что кристаллическая структура полученная экспериментаторами, предположительно, находилась в открытом состоянии. Было перепробовано множество попыток, прежде чем пришло понимание, что происходит "инактивация" фильтра это было нечто, что происходило спонтанно и препятствовало проводимости. В конце концов Халили-Арагхи поняла, что некоторые карбонатные группы в фильтре ответственны за непроводимость, и после многих попыток она ограничила их, после чего канал начал пропускать ионы. И, наконец, она и Шультен подтвердили, что проникновение ионов калия происходит по цепочному механизму похожему на перебор четок; в канале есть два иона калия, а затем входит третий, он остается там, но выталкивает нижний калий вниз через остальную часть канала. Хотя это подтверждение было опубликовано в Biophysical Letters, Халили-Арагхи говорит, что трудности моделирования проводимости были даже не самой сложной частью проекта.

Как-то она была на собрании Биофизического общества и завела разговор с ученым из лаборатории Бенуа Ру в Чикаго. Ру биофизик, который также изучал калиевые каналы. Халили-Арагхи узнал, что соавтор Ру, Владимир Яров, построил модель замкнутой структуры калиевого канала, закрытого напряжением. Это вскоре стало ценным для Халили-Арагхи, и вскоре они с Шультеном начали сотрудничать с лабораторией Ру. В то время как мембранные белки традиционно трудно кристаллизовать, еще труднее кристаллизовать ионный канал в состоянии, которое существует только при приложении определенного напряжения, чтобы держать его закрытым ибо это состояние должно происходить внутри мембраны. Но модель замкнутого состояния от Ярова не была похожа на первозданную кристаллографическую структуру. Халили-Арагхи пришлось потратить много времени на усовершенствование структуры, пока она не достигла стабильного закрытого состояния.

На этом работа не закончилась. "У нас было закрытое и открытое состояния, и мы должны были как-то проверить их, что закрытое состояние это действительно закрытое состояние", резюмирует она. Затем последовали серьезные вычисления расчет стробирующих зарядов. Когда канал переходит из открытого состояния в закрытое, некоторые заряженные остатки перемещаются, и эти движения могут быть количественно оценены в экспериментах. Хотя Халили-Арагхи могла точно воспроизвести эти экспериментальные измерения, сказать, что моделирование требовало огромного компьютерного времени, было бы преуменьшением. "В тот момент я думаю, что количество компьютерного времени, которое мы использовали для этого проекта, было, вероятно, сопоставимо суммарному затраченному времени на все другие проекты", говорит Халили-Арагхи. Но усилия того стоили, потому что они получили очень близкое совпадение с экспериментальными стробирующими зарядами. Халили-Арагхи провела еще много вычислений, чтобы проверить и улучшить закрытое состояние. В конце концов они обнаружили структуру, которая имела область сужения с отрицательно заряженными остатками. "Таким образом, все становилось на свои места, снова указывая на то, что наше закрытое состояние, вероятно, чем-то похоже на реальное закрытое состояние, обнаруженное в клеточной мембране нейрона", утверждает она. "Это, вероятно, самое близкое, что мы можем получить, я имею в виду, пока мы не увидим это с помощью кристаллографии или любого другого метода."

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

Тетрамерная структура калиевого каналаТетрамерная структура калиевого канала

Этот математически интересный способ управления нервными сигналами был открыт в 1950-х годах Ходжкином и Хаксли задолго до того, как стали известны какие-либо молекулярные подробности об ионных каналах в нейронной мембране. Ходжкин и Хаксли были удостоены в 1963 году Нобелевской премии по физиологии и медицине. Самое удивительное, что математическая форма уравнений Ходжкина-Хаксли напоминает недавно обнаруженную структуру ионных каналов, ясно указывая, что ионные каналы должны иметь четыре вентиля напряжения. До тех пор, пока Шультен не увлекся калиевыми каналами, оставалось загадкой, каким образом вентили напряжения улавливают слабые изменения напряжения, происходящие на нейронной мембране. Моделирование Халили-Арагхи и Шультена показало в случае калиевого канала, что изменение напряжения локализуется на удивительно коротком расстоянии в затворе, так что, хотя это изменение невелико, возникающая сила (которая зависит от времени изменения напряжения заряда иона и расстояния) достаточно высока, чтобы произвести сильный эффект. Физика включает в себя перераспределение высоко поляризуемых молекул воды и изменение обычной альфа-спирали белка в так называемую спираль 3-10. Эволюция открыла здесь электронное наноустройство, подобное транзистору, которое делает возможными нейронные сигналы и, следовательно, нейронные вычисления. В настоящее время группа Шультена расширяет исследования Халили-Арагхи на натриевый канал, который является более сложным по своему устройству. Шультен уверен, что вскоре он завершит мост между структурой и динамикой нейронных ионных каналов, с одной стороны, и математической феноменологией Ходжкина и Хаксли, с другой.

Грандиозный финал: серия взаимосвязанных процессов

Мы сотканы из ткани наших снов Шекспир, Буря

Две мечты Клауса Шультена определили его цели на протяжении последних сорока пяти лет научной жизни. В детстве он был очарован АТФ, и он никогда не забывал это раннее очарование универсальной валютой живой клетки. В том же духе он хотел изучить ключевой мембранный белок, реакционный центр, чтобы объяснить фотосинтез, и его возмутительный план построить свой собственный суперкомпьютер в 1980-х годах направил его на путь становления вычислительным биологом. Итак, спустя четыре десятилетия удовлетворен ли Шультен своим личным пониманием АТФ или фотосинтеза? Он подобрался очень близко. И чтобы представить кульминацию жизненного призвания, он счел необходимым осветить ряд процессов фотосинтеза, которые начинаются с поглощения солнечного света и заканчиваются выработкой АТФ. "И дело в том, что то, что я хотел сделать с помощью науки, я мог сделать только через фильм", резюмирует Шультен. Но это не голливудская постановка. Каждый отдельный атом правильно учтен, и все процессы были тщательно демистифицированы на протяжении многих лет многими, многими учеными, включая различных исследователей из группы Шультена. Никаких спецэффектов, только чистая наука. Этот ролик, предоставленный ниже, также затрагивает другую мечту Шультена, его желание описать, как несколько белков работают вместе в клетке, образуя сообщества биологической организации, что является сутью создания живой клетки. И этот фильм основан на мембранной системе, хроматофоре фиолетовых бактерий.

Хроматофор это автономная единица в фиолетовой бактерии, которая действует как биологическая машина, производя АТФ под действием солнечного света. На самом деле ее называют псевдо-органеллой, поскольку она встречается у прокариот. В принципе, он выполняет автономную функцию сам по себе, и хотя он действует только в контексте клетки, он имеет такие размеры, что теперь можно моделировать эту "органеллу" во всей ее полноте, и тем самым это приближает Шультена к цели изучения целой клетки. Хотя существует много различных форм хроматофоров, один из них, показанный в видео, имеет сферическую форму, около 60 нанометров в диаметре и содержит более 100 белков, встроенных в сферическую мембрану. Внутри хроматофора происходит ряд процессов между взаимодействующими белками. "Они работают, как швейцарские часы, где маленькие колесики и шестеренки приводят друг друга в движение", поясняет Шультен. "И именно этот часовой механизм процессов действительно заставляет систему делать то, что она делает." В то время как ученые в течение десятилетия рассматривали отдельные процессы, выполняемые отдельными белками, целью фильма Шультена было показать последовательность процессов, которые приводят к созданию АТФ из солнечного света, и проиллюстрировать взаимосвязь всех процессов.

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

Но как Шультен прошел путь от увлечения реакционным центром (одним из пяти типов белков, встроенных в хроматофор) до понимания всех типов белков, составляющих хроматофор? Как только структура реакционного центра была получена в 1985 году Гартмутом Мишелем, Шультен немедленно приступил к вычислениям этого мембранного белка, лежащего в основе фотосинтеза. И тогда он взялся за рискованный проект. Однажды вечером в Урбане в 1995 году на вечеринке Шультен и Мишель начали обсуждать смелый план по определению структуры светосборочного белка. В то время как реакционный центр имеет решающее значение для фотосинтеза, существует ряд белков, собирающих свет (показаны красным и зеленым на изображении), которые фактически направляют энергию в реакционный центр, где она затем преобразуется в перенос электронов. Белки собирающие свет, в основном действуют как вспомогательные антенны, чтобы поглощать больше доступного света для использования хроматофором. Пакт Шультена с Гартмутом Мишелем был на самом деле сумасшедшим, потому что Шультен не был экспертом в структурной биологии, а кристаллографические данные Гартмута упускали ключевую часть информации (фазовые углы). Но Шультен вместе с бесстрашным постдоком Сиче Ху фактически использовал вычисления, чтобы в конечном итоге выяснить структуру белка, называемого светосборочным комплексом. Таким образом, Шультен начал расширять свои знания о белковых компонентах хроматофора.

Моделирование LHC II (light harvesting complex) привело к тому, что Шультен называет своим самым большим успехом в области исследований мембранных белков, над которыми он работал с начала 1990-х годов. В тот момент Шультен снял фуражку структурного биолога и надел шляпу квантового физика. С выдающимся дуэтом аспирантов, Торстеном Ритцем и Анной Дамьянович, а также Сиче Ху, команда разгадала, как каротиноиды и хлорофиллы, составляющие светосборный комплекс II, функционируют как единое целое и используют квантовую физику для достижения эффективного поглощения фотонов и передачи энергии в реакционный центр. С одной стороны, Шультен был биофизиком, который занимался моделированием белков. "С другой стороны, я был физиком-теоретиком, который хорошо знал квантовую физику, и поэтому в тот момент, внезапно эта моя вторая половина мозга проснулась. И это было, я думаю, самое замечательное, что случилось со мной."

Таким образом, первоначально Шультен был очарован реакционным центром, а затем в середине 1990-х годов он начал работать над двумя типами вспомогательных белков, которые направляют энергию в реакционный центр. Кроме того, примерно в 2005 году Шультен начал амбициозный проект с постдоком Мелихом Сенером, чтобы собрать воедино все экспериментальные данные за последние десятилетия и фактически смоделировать полный хроматофор, вплоть до каждого последнего атома. Эта работа стала возможной в значительной степени благодаря сотруднику Шультена в Шеффилде, Нилу Хантеру. Было бы невозможно проиллюстрировать все взаимосвязанные процессы в хроматофоре до тех пор, пока не будет реализована его полная структура.

В хроматофоре есть пять ключевых белков-игроков. Первые два (светосборочные комплексы I и II) действуют как антенны, а третий белок, реакционный центр, преобразует собранную энергию в перенос электронов. Таким образом, остаются два типа белков, необходимых для окончательного превращения солнечного света в АТФ. Эти два комплекса называются bc1-комплексом и АТФ-синтазой. После того, как свет попадает на антенны хроматофора, он поступает в реакционный центр, который использует эту энергию, чтобы протолкнуть электрон вниз по цепи. Затем специальная молекула (хинон) активируется двумя электронами и двумя протонами (становится хинолом) и проплывает через маслянистую мембрану к bc1-комплексу. Затем bc1 перекачивает протоны хинола с внешней стороны сферического хроматофора внутрь. И тогда хорошо известная АТФ-синтаза использует протонный градиент, чтобы стимулировать производство АТФ. (Для полноты картины, здесь задействован шестой белок цитохром, который переносит электроны от bc1 обратно к реакционному центру.)

Группа Шультена должна была использовать инструменты квантовой химии для описания сложных реакций переноса электронов и протонов, которые происходят в bc1-комплексе, усилия, возглавляемые совсем недавно аспиранткой Анджелой Барраган. Она провела очень тщательные крупномасштабные квантово-химические расчеты, чтобы точно определить, что на самом деле происходит в процессе создания протонного градиента.

Как и следовало ожидать, работа группы Шультена над АТФ-синтазой дорога его сердцу. Ему потребовалось много лет, чтобы достичь подходящего понимания этого мембранного белка. Поначалу все шло медленно, поскольку первоначальные структуры не были достаточно полными, чтобы Шультен мог полностью погрузиться в объяснение его внутренней работы. "Поэтому тогда я был очень осторожен и опубликовал только то немногое, что понял", рассказывает Шультен о своей многолетней работе над АТФ-синтазой. "И сегодня я, наконец, думаю, что у меня все получилось." Причина, по которой Шультен считает, что он, наконец, может объяснить АТФ-синтазу и ее функцию, затрагивает тему его фильма о хроматофоре в целом: совместная работа всех взаимосвязанных частей. АТФ-синтаза содержит шесть белков. Только сегодня, благодаря работе, недавно выполненной фантастическим постдоком Абхи Сингхарой, Шультен может описать сложность системы, а именно, как шесть белков работают вместе, чтобы преобразовать протонный градиент в крутящий момент, а затем в настоящую АТФ. Хотя он все еще проверяет некоторые шаги, используемые АТФ-синтазой, он находится в пределах досягаемости своих детских мечт об АТФ.

YouTube аналог

Фильм "Хроматофор", таким образом, заключает в себе множество тем, которые извивались в мозгу Шультена в течение последних пятидесяти лет. Он с детства хотел исследовать молекулы жизни, особенно важнейшую АТФ. Поэтому он принял решение изучать молекулярную физику в аспирантуре Гарварда дисциплину, которую все его друзья считали скучной и недостойной внимания. Во время своей первой работы он был знаком с областью фотосинтеза, но мог проводить только ограниченные исследования на эту тему, поскольку не было намеков, что важнейшие мембранные белки, лежащие в основе фотосинтеза, будут когда-либо оцифрованны. В Мюнхене в 1980-х годах Шультен имел место в первом ряду для получения структуры первого мембранного белка, реакционного центра. Но реакционный центр имел и другой глубокий эффект, поскольку он подтолкнул Шультена к использованию компьютера в качестве инструмента, который, если его отточить, может стать бесценным для биологии. И, наконец, он всегда настаивал на моделировании все больших и больших структур с явной целью описания биологической организации, что переводится в рассмотрение того, как взаимодействуют кластеры из сотен или даже тысяч белков. И у него были программные средства, которые он разработал для вычислительной биологии, VMD и NAMD, и с помощью этих самых инструментов он был способен сделать фильм про хроматофор, который требовал моделирования и анализа огромного числа атомов (около 100 миллионов). Имея все это в виду, Шультен утверждает, что: "Идея о фильме напрашивалась сама собой."

Фильм "Хроматофор" является подходящим завершением этой истории о мембранных белках. С богатым набором мембранных белков, которые Шультен и его группа уже демистифицировали, удивительно, что Шультену пришлось ждать почти тридцать лет, пока он смог смоделировать белок и мембрану вместе. Правда, его путь не всегда был легким, но Шультен смотрел вперед и не зацикливался на невзгодах. Его успехи в этой области идут рука об руку с достижениями в области вычислительных инструментов, которые разрабатывает его группа. На самом деле, с помощью компьютера Шультен неплохо зарабатывал себе на жизнь, хотя получение структур мембранных белков с помощью кристаллографии все еще происходит довольно медленно. Хотя роль вычислительного биолога иногда недооценивается, Шультен любит свое жизненное призвание. "Итак, теперь я разбираюсь в живых системах с помощью своего математического ума, и я очень, очень доволен этим".


Статья была написана в 2014 г. Через год хроматофор (100 миллионов атомов) промоделировали на суперкомпьютере "Титан". В октябре 2016-го Клаусс Шультен умер в возрасте 69 лет, но его дело продолжается его учениками из Theoretical and Computational Biophysics Group.

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

Подробнее..

Перенос молекулярной динамики на CUDA. Часть II Суммирование по Эвальду

07.07.2020 20:09:59 | Автор: admin
В предыдущей статье мы обсудили основу метода молекулярной динамики, в том числе вычисление энергии и сил взаимодействия между частицами с заданными парными потенциалами. А что, если частицы обладают некоторым электрическим зарядом? Например, в том случае, если мы моделируем кристалл поваренной соли, состоящий из ионов Na+ и Cl-. Или водный раствор, содержащий те или иные ионы. В этом случае, кроме парных потенциалов типа Леннарда-Джонса между ионами действуют силы электростатического взаимодействия, т.е. закон Кулона. Энергия такого взаимодействия для пары частиц i-j равна:

$E=C\frac{q_iq_j}{r_{ij}},$


где q заряд частицы, rij расстояние между частицами, С некоторая постоянная, зависящая от выбора единиц измерения. В системе СИ это $\frac{1}{4\pi\epsilon_0}$, в СГС 1, в моей программе (где энергия выражена в электронвольтах, расстояние в ангстремах, а заряд в элементарных зарядах) C примерно равно 14.3996.
image
Ну и что, скажете вы? Просто добавим соответствующее слагаемое в парный потенциал и готово. Однако, чаще всего в МД моделировании используют периодические граничные условия, т.е. моделируемая система со всех сторон окружена бесконечным количеством её виртуальных копий. В этом случае каждый виртуальный образ нашей системы будет взаимодействовать со всеми заряженными частицами внутри системы по закону Кулона. А поскольку Кулоновское взаимодействие убывает с расстоянием очень слабо (как 1/r), то отмахнуться от него так просто нельзя, сказав, что с такого-то расстояния мы его не вычисляем. Ряд вида 1/x расходится, т.е. его сумма, в принципе, может расти до бесконечности. И что же теперь, миску супа не солить? Убьёт электричеством?

можно не только солить суп, но и посчитать энергию Кулоновского взаимодействия в периодических граничных условиях. Такой метод был предложен Эвальдом ещё в 1921 году для расчета энергии ионного кристалла (можно также посмотреть в википедии). Суть метода заключатся в экранировании точечных зарядов и последующим вычетом функции экранирования. При этом часть электростатического взаимодействия сводится к короткойдействующему и его можно просто обрезать стандартным способом. Оставшаяся дальнодействующая часть эффективно суммируется в пространстве Фурье. Опуская вывод, который можно посмотреть в статье Блинова или в той же книге Френкеля и Смита сразу запишу решение, называемое суммой Эвальда:

$E_{real} = С\sum\limits^{N-1}_i{\sum\limits^N_{j=i+1}{ \frac{q_iq_j}{r_{ij}}\operatorname{erfc}(\alpha r_{ij}) }}, $


$E_{rec} = С\frac{2\pi}{V}\sum_{\vec{k}\ne0}(\exp(-k^2/4\alpha^2)/k^2 \sum_jq_i\exp(i\vec{k}\vec{r}_j)),$


$E_{const} = -\frac{C}{V\alpha} \sum_iq_i^2, $


где параметр, регулирующий соотношение вычислений в прямом и обратном пространствах, k вектора в обратном пространстве по которым идёт суммирование, V объём системы (в прямом пространстве). Первая часть (Ereal) является короткодействующей и вычисляется в том же цикле, что и другие парные потенциалы, смотри функцию real_ewald в предыдущей статье. Последний вклад (Eсonst) является поправкой на самовзаимодействие и часто называется постоянной частью, поскольку не зависит от координат частиц. Её вычисление тривиально, поэтому мы остановимся только на второй части Эвальдовой суммы (Erec), суммировании в обратном пространстве. Естественно, во времена вывода Эвальда молекулярной динамики не было, кто впервые использовал этот метод в МД найти мне не удалось. Сейчас любая книга по МД содержит его изложение как некий золотой стандарт. К книге Аллена даже прилагается пример кода на фортране. К счастью, у меня остался код, написанный когда-то на С для последовательной версии, осталось только его распараллелить (я позволил себе опустить некоторые объявления переменных и другие несущественные детали):
void ewald_rec(){    int mmin = 0;    int nmin = 1;    // массивы где хранятся iexp(x[i] * kx[l]),    double** elc;    double** els;    //... iexp(y[i] * ky[m]) и    double** emc;    double** ems;    //... iexp(z[i] * kz[n]),    double** enc;    double** ens;    // временные массивы для произведений iexp(x*kx)*iexp(y*ky)    double* lmc;    double* lms;    // и для q[i] * iexp(x*kx)*iexp(y*ky)*iexp(z*kz)    double* ckc;    double* cks;    // ПРЕДВАРИТЕЛЬНОЕ ЗАПОЛНЕНИЕ МАССИВОВ    eng = 0.0;    for (i = 0; i < Nat; i++)   // цикл по атомам    {        // emc/s[i][0] и enc/s[i][0] уже заполнены на этапе инициализации        // в массив elc/s нужно обновить, см. далее        elc[i][0] = 1.0;        els[i][0] = 0.0;        // iexp(kr)        sincos(twopi * xs[i] * ra, els[i][1], elc[i][1]);        sincos(twopi * ys[i] * rb, ems[i][1], emc[i][1]);        sincos(twopi * zs[i] * rc, ens[i][1], enc[i][1]);    }    // заполняем следующие элементы массива emc/s[i][l] = iexp(y[i]*ky[l]) итеративно, используя комплексное умножение    for (l = 2; l < ky; l++)        for (i = 0; i < Nat; i++)        {            emc[i][l] = emc[i][l - 1] * emc[i][1] - ems[i][l - 1] * ems[i][1];            ems[i][l] = ems[i][l - 1] * emc[i][1] + emc[i][l - 1] * ems[i][1];        }    // заполняем следующие элементы массива enc/s[i][l] = iexp(z[i]*kz[l]) итеративно, используя комплексное умножение    for (l = 2; l < kz; l++)        for (i = 0; i < Nat; i++)        {            enc[i][l] = enc[i][l - 1] * enc[i][1] - ens[i][l - 1] * ens[i][1];            ens[i][l] = ens[i][l - 1] * enc[i][1] + enc[i][l - 1] * ens[i][1];        }    // ГЛАВНЙ ЦИКЛ ПО ВСЕМ K-ВЕКТОРАМ:    for (l = 0; l < kx; l++)    {        rkx = l * twopi * ra;         // записываем exp(ikx[l]) в ikx[0] для сохранения памяти        if (l == 1)            for (i = 0; i < Nat; i++)            {                elc[i][0] = elc[i][1];                els[i][0] = els[i][1];            }        else if (l > 1)            for (i = 0; i < Nat; i++)            {                // iexp(kx[0]) = iexp(kx[0]) * iexp(kx[1])                x = elc[i][0];                elc[i][0] = x * elc[i][1] - els[i][0] * els[i][1];                els[i][0] = els[i][0] * elc[i][1] + x * els[i][1];            }        for (m = mmin; m < ky; m++)        {            rky = m * twopi * rb;            // заполняем временный массив произведением iexp(kx*x[i]) * iexp(ky*y[i])            if (m >= 0)                for (i = 0; i < Nat; i++)                {                    lmc[i] = elc[i][0] * emc[i][m] - els[i][0] * ems[i][m];                    lms[i] = els[i][0] * emc[i][m] + ems[i][m] * elc[i][0];                }            else // для отрицательных значений m используем комплексное сопряжение:                for (i = 0; i < Nat; i++)                {                    lmc[i] = elc[i][0] * emc[i][-m] + els[i][0] * ems[i][-m];                    lms[i] = els[i][0] * emc[i][-m] - ems[i][-m] * elc[i][0];                }            for (n = nmin; n < kz; n++)            {                rkz = n * twopi * rc;                rk2 = rkx * rkx + rky * rky + rkz * rkz;                if (rk2 < rkcut2) // используем радиус обрезания                {                    // вычисляем сумму (q[i]*iexp(kr[k]*r[i])) - зарядовую плотность                    sumC = 0; sumS = 0;                    if (n >= 0)                        for (i = 0; i < Nat; i++)                        {                            //считываем заряд частицы                            ch = charges[types[i]].charge;                            ckc[i] = ch * (lmc[i] * enc[i][n] - lms[i] * ens[i][n]);                            cks[i] = ch * (lms[i] * enc[i][n] + lmc[i] * ens[i][n]);                            sumC += ckc[i];                            sumS += cks[i];                        }                    else // для отрицательных индексов используем комплексное сопряжение:                        for (i = 0; i < Nat; i++)                        {                            //считываем заряд частицы                            ch = charges[types[i]].charge;                            ckc[i] = ch * (lmc[i] * enc[i][-n] + lms[i] * ens[i][-n]);                            cks[i] = ch * (lms[i] * enc[i][-n] - lmc[i] * ens[i][-n]);                            sumC += ckc[i];                            sumS += cks[i];                        }                    //наконец вычисляем энергию и силы                    akk = exp(rk2 * elec->mr4a2) / rk2;                    eng += akk * (sumC * sumC + sumS * sumS);                    for (i = 0; i < Nat; i++)                    {                        x = akk * (cks[i] * sumC - ckc[i] * sumS) * C * twopi * 2 * rvol;                        fxs[i] += rkx * x;                        fys[i] += rky * x;                        fzs[i] += rkz * x;                    }                }            } // end n-loop (over kz-vectors)            nmin = 1 - kz;        } // end m-loop (over ky-vectors)        mmin = 1 - ky;    }  // end l-loop (over kx-vectors)   engElec2 = eng * С * twopi * rvol;}

Пара пояснений к коду: функция считает комплексную экспоненту (в комментариях к коду она обозначена iexp, чтобы убрать мнимую единицу из скобок) от векторного произведения k-вектора на радиус-вектор частицы для всех k-векторов и для всех частиц. Эта экспонента домножается на заряд частицы. Далее вычисляется сумма таких произведений по всем частицам (внутренняя сумма в формуле для Erec), у Френкеля она называется зарядовой плотностью, а у Блинова структурным фактором. Ну а далее, на основании этих структурных факторов вычисляется энергия и силы, действующие на частицы. Компоненты k-векторов (2*l/a, 2*m/b, 2*n/c) характеризуются тройкой целых чисел l, m и n, которые и пробегаются в циклах до заданных пользователем пределов. Параметры a, b и c это размеры моделируемой системы в измерениях x, y и z соответственно (вывод верен для системы с геометрией прямоугольного параллелепипеда). В коде 1/a, 1/b и 1/с соответствуют переменным ra, rb и rc. Массивы под каждую величину представлены в двух экземплярах: под действительную и мнимую части. Каждый следующий k-вектор в одном измерении получается итеративно из предыдущего путем комплексного умножения предыдущего на единичный, чтобы каждый раз не считать синус с косинусом. Массивы emc/s и enc/s заполняются для всех m и n, соответственно, а массив elc/s значение для каждого l>1 помещает в нулевой индекс по l в целях экономии памяти.
В целях распараллеливания выгодно вывернуть порядок циклов так, чтобы внешний цикл пробегался по частицам. И тут мы видим проблему распараллелить эту функцию можно только до вычисления суммы по всем частицам (зарядовой плотности). Дальнейшие вычисления опираются на эту сумму, а она будет рассчитана только когда все потоки закончат работу, поэтому придётся разбить эту функцию на две. Первая вычисляет считает зарядовую плотность, а вторая энергию и силы. Замечу, что во второй функции снова потребуется величина qiiexp(kr) для каждой частицы и для каждого k-вектора, вычисленная на предыдущем этапе. И тут есть два подхода: либо пересчитать её заново, либо запомнить. Первый вариант требует больше времени, второй больше памяти (количество частиц * количество k-векторов * sizeof(float2)). Я остановился на втором варианте:
__global__ void recip_ewald(int atPerBlock, int atPerThread, cudaMD* md)// calculate reciprocal part of Ewald summ// the first part : summ (qiexp(kr)) evaluation{    int i;      // for atom loop    int ik;     // index of k-vector    int l, m, n;    int mmin = 0;    int nmin = 1;    float tmp, ch;    float rkx, rky, rkz, rk2;   // component of rk-vectors    int nkx = md->nk.x;    int nky = md->nk.y;    int nkz = md->nk.z;        // arrays for keeping iexp(k*r) Re and Im part    float2 el[2];    float2 em[NKVEC_MX];    float2 en[NKVEC_MX];    float2 sums[NTOTKVEC];          // summ (q iexp (k*r)) for each k    extern __shared__ float2 sh_sums[];     // the same in shared memory    float2 lm;     // temp var for keeping el*em    float2 ck;     // temp var for keeping q * el * em * en (q iexp (kr))    // invert length of box cell    float ra = md->revLeng.x;    float rb = md->revLeng.y;    float rc = md->revLeng.z;    if (threadIdx.x == 0)        for (i = 0; i < md->nKvec; i++)            sh_sums[i] = make_float2(0.0f, 0.0f);    __syncthreads();    for (i = 0; i < md->nKvec; i++)        sums[i] = make_float2(0.0f, 0.0f);    int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    ik = 0;    for (i = id0; i < N; i++)    {        // save charge        ch = md->specs[md->types[i]].charge;        el[0] = make_float2(1.0f, 0.0f);    // .x - real part (or cos) .y - imagine part (or sin)                em[0] = make_float2(1.0f, 0.0f);        en[0] = make_float2(1.0f, 0.0f);        // iexp (ikr)        sincos(d_2pi * md->xyz[i].x * ra, &(el[1].y), &(el[1].x));        sincos(d_2pi * md->xyz[i].y * rb, &(em[1].y), &(em[1].x));        sincos(d_2pi * md->xyz[i].z * rc, &(en[1].y), &(en[1].x));        // fil exp(iky) array by complex multiplication        for (l = 2; l < nky; l++)        {             em[l].x = em[l - 1].x * em[1].x - em[l - 1].y * em[1].y;             em[l].y = em[l - 1].y * em[1].x + em[l - 1].x * em[1].y;        }        // fil exp(ikz) array by complex multiplication        for (l = 2; l < nkz; l++)        {             en[l].x = en[l - 1].x * en[1].x - en[l - 1].y * en[1].y;             en[l].y = en[l - 1].y * en[1].x + en[l - 1].x * en[1].y;        }        // MAIN LOOP OVER K-VECTORS:        for (l = 0; l < nkx; l++)        {            rkx = l * d_2pi * ra;                         // move exp(ikx[l]) to ikx[0] for memory saving (ikx[i>1] are not used)            if (l == 1)                el[0] = el[1];            else if (l > 1)                {                    // exp(ikx[0]) = exp(ikx[0]) * exp(ikx[1])                    tmp = el[0].x;                    el[0].x = tmp * el[1].x - el[0].y * el[1].y;                    el[0].y = el[0].y * el[1].x + tmp * el[1].y;                }            //ky - loop:            for (m = mmin; m < nky; m++)            {                rky = m * d_2pi * rb;                //set temporary variable lm = e^ikx * e^iky                if (m >= 0)                {                        lm.x = el[0].x * em[m].x - el[0].y * em[m].y;                         lm.y = el[0].y * em[m].x + em[m].y * el[0].x;                }                else // for negative ky give complex adjustment to positive ky:                {                        lm.x = el[0].x * em[-m].x + el[0].y * em[-m].y;                        lm.y = el[0].y * em[-m].x - em[-m].x * el[0].x;                }                //kz - loop:                for (n = nmin; n < nkz; n++)                {                    rkz = n * d_2pi * rc;                    rk2 = rkx * rkx + rky * rky + rkz * rkz;                    if (rk2 < md->rKcut2) // cutoff                    {                        // calculate summ[q iexp(kr)]   (local part)                        if (n >= 0)                         {                                ck.x = ch * (lm.x * en[n].x - lm.y * en[n].y);                                ck.y = ch * (lm.y * en[n].x + lm.x * en[n].y);                         }                        else // for negative kz give complex adjustment to positive kz:                         {                                ck.x = ch * (lm.x * en[-n].x + lm.y * en[-n].y);                                ck.y = ch * (lm.y * en[-n].x - lm.x * en[-n].y);                        }                        sums[ik].x += ck.x;                        sums[ik].y += ck.y;                                                // save qiexp(kr) for each k for each atom:                        md->qiexp[i][ik] = ck;                        ik++;                    }                } // end n-loop (over kz-vectors)                nmin = 1 - nkz;            } // end m-loop (over ky-vectors)            mmin = 1 - nky;        }  // end l-loop (over kx-vectors)    } // end loop by atoms    // save sum into shared memory    for (i = 0; i < md->nKvec; i++)    {        atomicAdd(&(sh_sums[i].x), sums[i].x);        atomicAdd(&(sh_sums[i].y), sums[i].y);    }    __syncthreads();    //...and to global    int step = ceil((double)md->nKvec / (double)blockDim.x);    id0 = threadIdx.x * step;    N = min(id0 + step, md->nKvec);    for (i = id0; i < N; i++)    {        atomicAdd(&(md->qDens[i].x), sh_sums[i].x);        atomicAdd(&(md->qDens[i].y), sh_sums[i].y);    }}// end 'ewald_rec' function__global__ void ewald_force(int atPerBlock, int atPerThread, cudaMD* md)// calculate reciprocal part of Ewald summ// the second part : enegy and forces{    int i;      // for atom loop    int ik;     // index of k-vector    float tmp;    // accumulator for force components    float3 force;    // constant factors for energy and force    float eScale = md->ewEscale;    float fScale = md->ewFscale;    int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    for (i = id0; i < N; i++)    {        force = make_float3(0.0f, 0.0f, 0.0f);        // summ by k-vectors        for (ik = 0; ik < md->nKvec; ik++)        {            tmp = fScale * md->exprk2[ik] * (md->qiexp[i][ik].y * md->qDens[ik].x - md->qiexp[i][ik].x * md->qDens[ik].y);            force.x += tmp * md->rk[ik].x;            force.y += tmp * md->rk[ik].y;            force.z += tmp * md->rk[ik].z;        }        md->frs[i].x += force.x;        md->frs[i].y += force.y;        md->frs[i].z += force.z;    } // end loop by atoms    // one block calculate energy    if (blockIdx.x == 0)        if (threadIdx.x == 0)        {            for (ik = 0; ik < md->nKvec; ik++)                md->engCoul2 += eScale * md->exprk2[ik] * (md->qDens[ik].x * md->qDens[ik].x + md->qDens[ik].y * md->qDens[ik].y);        }}// end 'ewald_force' function

Надеюсь, вы мне простите, что я оставил комментарии на английском, код практически повторяет последовательную версию. Код даже стал читабельнее за счет того, что массивы потеряли одно измерение: elc/s[i][l], emc/s[i][m] и enc/s[i][n] превратились в одномерные el, em и en, массивы lmc/s и ckc/s в переменные lm и ck (пропала мерность по частицам, поскольку отпала необходимость хранить это для каждой частицы, промежуточный результат накапливается в shared memory). К сожалению, тут же возникла и проблема: массивы em и en пришлось задать статическими, чтобы не использовать глобальную память и не выделять память динамически каждый раз. Количество элементов в них определяется директивой NKVEC_MX (максимальное количество k-векторов по одному измерению) на этапе компиляции, а runtime используются только первые nky/z элементов. Кроме того, появился сквозной индекс по всем k-векторам и аналогичная директива, ограничивающая общее количество этих векторов NTOTKVEC. Проблема возникнет, если пользователю понадобится больше k-векторов, чем определено директивами. Для вычисления энергии предусмотрен блок с нулевым индексом, поскольку неважно какой именно блок выполнит этот расчет и в каком порядке. Тут может быть надо было использовать готовые функции FFT, раз уж метод основан на нём, но я так и не сообразил, как это сделать.
Ну а теперь попробуем что-нибудь посчитать, да тот же содиум хлорайд. Возьмём 2 тысячи ионов натрия и столько же хлора. Заряды зададим целыми, а парные потенциалы возьмём, например, из этой работы. Стартовую конфигурацию зададим случайно и слегка перемешаем её, рисунок 2а. Объём системы выберем так, чтобы он соответствовал плотности поваренной соли при комнатной температуре (2,165 г/см3). Запустим все это на небольшое время (10000 шагов по 5 фемтосекунд) с наивным учетом электростатики по закону Кулона и используя суммирование по Эвальду. Результирующие конфигурации приведены на рисунках 2б и 2в, соответственно. Вроде бы в случае с Эвальдом система чуть больше упорядочилась чем без него. Важно также, что флуктуации полной энергии с применением суммирования существенно уменьшились.

Рисунок 2. Начальная конфигурация системы NaCl (a) и после 10000 шагов интегрирования: наивным способом (б) и со схемой Эвальда (в).

Вместо заключения


Замечу, что структура, получаемая на рисунке, не соответствует кристаллической решетки NaCl, а скорее решетке ZnS, но это уже претензия к парным потенциалам. Учет же электростатики очень важен для молекулярно-динамического моделирования. Считается, что именно электростатическое взаимодействие ответственно за образование кристаллических решёток, поскольку действует на больших расстояниях. Правда с этой позиции сложно объяснить как при охлаждении кристаллизуются такие вещества как аргон.
Кроме упомянутого метода Эвальда, есть ещё и другие способы учета электростатики, смотрите, например, этот обзор.
Подробнее..

Перенос молекулярной динамики на CUDA. Часть III Внутримолекулярное взаимодействие

21.08.2020 16:19:43 | Автор: admin
До этого мы рассматривали молекулярную динамику, где законы взаимодействия между частицами зависели исключительно от типа частиц или от их заряда. Для веществ молекулярной природы взаимодействие между частицами (атомами) сильно зависит от того, принадлежат ли атомы одной молекуле или нет (точнее, связаны ли они химической связью). Например, вода:
image
очевидно, что водород с кислородом внутри одной молекулы взаимодействуют совсем по-другому, нежели тот же кислород с водородом соседней молекулы. Таким образом, разделяют ВНУТРИмолекулярное (intramolecular) и МЕЖмолекулярное (intermolecular) взаимодействие. Межмолекулярное взаимодействие можно задать короткодействующими и Кулоновскими парными потенциалами, о которых речь шла в предыдущих статьях. Здесь же сконцентрируемся на внутримолекулярном.

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

$U=\frac12k(r-r_0)^2,$


где r расстояние между частицами, k константа жесткости связи, а r0 равновесная длина связи; и потенциал Морзе:

$U=D(1-\exp(-\alpha(r-r_0)))^2,$


где D глубина потенциальной ямы, параметр характеризует ширину потенциальной ямы.
Следующий тип внутримолекулярных взаимодействий это валентные углы. Обратимся снова к заглавной картинке. Почему молекула изображена уголком, ведь силы электростатики должны были обеспечить максимальное расстояние между ионами водорода, которое соответствует углу H-O-H равному 180? Дело в том, что на рисунке нарисовано не все. Из школьного курса химии, можно вспомнить, что у кислорода есть ещё 2 неподелённые электронные пары, взаимодействие с которыми и искажает угол:
image
В классической молекулярной динамике обычно не вводят такие объекты как электроны или электронные облака, поэтому для имитации правильных углов используют потенциал валентного угла, т.е. функциональную зависимость потенциальной энергии от координат 3х частиц. Одним из наиболее удобных таких потенциалов является гармонический косинус:

$U=\frac12k(\theta-\theta_0)^2,$


где угол образованной тройкой частиц, k константа жесткости, а 0 равновесный угол.
Существуют внутримолекулярные потенциалы и более высокого порядка, например, торсионные углы, но они ещё более искусственные, чем валентные углы.
Добавить взаимодействия между частицами с заранее известными индексами дело тривиальное. Для связей храним массив, содержащий индексы связанных частиц и тип взаимодействия. Каждому потоку выдаем свой диапазон связей на обработку и готово. Аналогично с валентными углами. Поэтому мы сразу себе усложним жизнь задачу: добавим возможность создавать/удалять химические связи и валентные углы runtime. Это сразу выводит нас из плоскости классической молекулярной динамики и открывает новый горизонт возможностей. Иначе можно было просто скачать что-нибудь из существующих пакетов, например LAMMPS, DL_POLY или GROMACS , тем более что они распространяются бесплатно.
Ну а теперь немного кода. Добавим соответствующие поля в основную структуру:
    //bonds:    int nBond;      // number of bonds    int mxBond;          // maximal number of bonds    int4* bonds;    // array of bonds     int* nbonds;    // count of bond for a given atom    int* neighToBind;   // a neighbor of a given atom for binding    int* canBind;       // flags that atom[iat] can be bind    int* r2Min;         // distances for the nearest neighbor (used for binding)    int* parents;       // indexes of one of the atom bonded with a given    cudaBond* bondTypes;     int** def_bonds;    // array[nSpec][nSpec] of default bond types    int** bindBonds;    // array[nSpec][nSpec] bond types created by binding    float** bindR2;        // square of binding distance [nSpec][nSpec]    //angles:    int nAngle;    // number of angles    int mxAngle;    int4* angles;   // array of angles      int* nangles;        // number of angles for given atom    int* oldTypes;          cudaAngle* angleTypes;     int* specAngles;    // [nSp] angle type formed by given species


Количество связей и углов переменное, но практически всегда можно оценить максимально возможное и выделить память сразу под максимум, чтобы не перевыделять память, поля nBond и mxBond, соответственно, означают текущее число связей и максимальное. Массив bonds будет содержать индексы связываемых атомов, тип связи и время создания связи (если нас вдруг заинтересует такая статистика, как среднее время жизни связи). Массив angles будет хранить индексы тройки атомов, образующих валентный угол, и тип валентного угла. Массивы bondTypes и angleTypes будут содержать характеристики возможных потенциалов валентных связей и углов. Вот их структуры:
struct cudaBond{    int type;  // potential type    int spec1, spec2; // type of atoms that connected by this bond type    int new_type[2];      // bond type after mutation    int new_spec1[2], new_spec2[2];    int mxEx, mnEx;     // flags: maximum or minimum of bond length exists    float p0, p1, p2, p3, p4;    // potential parameters    float r2min, r2max;          // square of minimal and maximal bond length    float (*force_eng)(float r2, float r, float &eng, cudaBond *bond); // return energy     int count;      // quantity of such bonds    float rSumm;        // summ of lentghs (for mean length calculation)    int rCount;          // number of measured lengths (for mean length calculation)    int ltSumm, ltCount;    // for calculation of lifetime};struct cudaAngle{    int type; // potential type    float p0, p1, p2;    // potential parameters    void (*force_eng)(int4* angle, cudaAngle* type, cudaMD* md, float& eng);};

Поле type определяет функциональный форму потенциала, new_type, new_spec1 и new_spec это индексы типа связи и типов связываемых атомов, после того как связь изменится (оборвется или превратится в связь другого типа). Эти поля представлены в виде массивов с двумя элементами. Первый соответствует ситуации, когда длина станет короче r2min1/2, второй когда превысит r2max1/2. Наиболее сложная часть алгоритма это применение свойств всех связей, с учетом возможности их обрыва и превращения, а также того, что другие потоки могли оборвать соседние связи, что привело к изменению типа связанных атомов. Поясню на примере той же воды. Изначально молекула электронейтральна, химические связи образованы электронами общими для водорода и кислорода. Грубо говоря можно сказать, что заряды на атомах водорода и кислорода нулевые (на самом деле электронная плотность смещена к кислороду, поэтому на водороде небольшой плюс, +, а на кислороде 2-). Если мы будем рвать связь, кислород окончательно заберет себе электрон, а водород его отдаст. Получатся частицы H+ и O-. Итого получается 5 типов частиц, условно обозначим их: H, H+, O, O-, O2-. Последняя образуется, если мы оторвем оба водорода от молекулы воды. Соответственно, реакции:
H2O -> H+ + OH-
и
OH- -> H+ + O2-.
Знатоки химии меня поправят, что при стандартных условиях для воды и первая то стадия распада практически не реализуется (в равновесии, только одна молекула из 107 диссоциирована на ионы, да и то не совсем такие, как написано). Но для описания алгоритмов, такие схемы будут наглядными. Допустим, один поток обрабатывает одну связь в молекуле воды, а другой поток вторую связь этой же молекулы. И так получилось, что обе связи нужно разорвать. Тогда один поток должен превратить атомы в H+ и O-, а второй в H+ и O2-. Но если потоки это делают одновременно, на момент начала процедуры кислород находится в состоянии O и оба потока его превращают в O-, что неверно. Вот такие ситуации нам необходимо как-то предотвращать. Блок-схема функции, обрабатывающей химическую связь:

Проверяем соответствуют ли текущие типы атомов типу связи, если нет, то берем из таблицы дефолтных типов (она должна быть заранее составлена), далее определяем квадрат расстояния между атомами (r2) и, если связь подразумевает максимальную или минимальную длину, проверяем, не вышли ли мы за эти границы. Если вышли, то нам необходимо поменять тип связи или удалить её и в обоих случаях поменять типы атомов. Для этого будет использована функция atomicCAS сравниваем текущий тип атома с тем, который должен быть и в этом случае заменяем на новый. Если тип атома уже был изменен другим потоком, возвращаемся в начало, к переопределению типа связи. Наихудший вариант, если тип первого атома мы успели поменять, а второго нет. Откатываться назад уже поздно, ведь после того, как мы поменяли первый атом, другие потоки уже могли что-то с ним сделать. Какой же выход? Предлагаю притвориться, что мы рвём/изменяем связь другого типа, а не того, за которую взялись вначале. Находим, какой тип связи должен был быть между начальным первым атомом и измененным вторым и обрабатываем её по тем же правилам, что и изначально ожидаемую. Если и в этом случае тип атома снова успел поменяться снова используем эту же схему. Тут неявно подразумевается, что новый тип связи обладает теми же свойствами рвется при той же длине и т.д., а частицы, образуемые в ходе разрыва такие, какие нужно. Поскольку эту информацию задет пользователь, мы как бы перекладываем ответственность с нашей программы на него, он должен задать все корректно. Код:
__global__ void apply_bonds(int iStep, int bndPerBlock, int bndPerThread, cudaMD* md){    int def;    int id1, id2;       // atom indexes    int old, old_spec2, spec1, spec2, new_spec1, new_spec2;     // atom types    int new_bond_type;        int save_lt, need_r, loop;    // flags to save lifetime, to need to calculate r^2 and to be in while loop    int mnmx;   // flag minimum or maximum    int action; // flag: 0 - do nothing, 1 - delete bond, 2 - transform bond    cudaBond *old_bnd, *cur_bnd;// old bond type, current bond type    float dx, dy, dz, r2, r;    float f, eng = 0.0f;    __shared__ float shEng;#ifdef DEBUG_MODE    int cnt;    // count of change spec2 loops#endif    if (threadIdx.x == 0)    {        shEng = 0.0f;    }    __syncthreads();    int id0 = blockIdx.x * bndPerBlock + threadIdx.x * bndPerThread;    int N = min(id0 + bndPerThread, md->nBond);    int iBnd;    for (iBnd = id0; iBnd < N; iBnd++)      if (md->bonds[iBnd].z)  // the bond is not broken      {          // atom indexes          id1 = md->bonds[iBnd].x;          id2 = md->bonds[iBnd].y;          // atom types          spec1 = md->types[id1];          spec2 = md->types[id2];          old_bnd = &(md->bondTypes[md->bonds[iBnd].z]);          cur_bnd = old_bnd;          save_lt = 0;          need_r = 1;          loop = 1;#ifdef DEBUG_MODE          cnt = 0;#endif                    if ((cur_bnd->spec1 == spec1)&&(cur_bnd->spec2 == spec2))          {              //ok          }          else              if ((cur_bnd->spec1 == spec2) && (cur_bnd->spec2 == spec1))              {                  invert_bond(id1, id2, spec1, spec2, &(md->bonds[iBnd]));                  //... then ok              }              else // atom types do not correspond to bond types              {                  save_lt = 1;              }          // end initial stage          while (loop)          {             if (save_lt)                    {                  def = md->def_bonds[spec1][spec2];                  if (def == 0)     // these atom types do not form a bond                  {#ifdef DEBUG_MODE                      printf("probably, something goes wrong\n");#endif                      action = 1;   // delete                      break;                  }                  else                  {                      //! change bond type and go on                      if (def < 0)                        {                          invert_bond(id1, id2, spec1, spec2, &(md->bonds[iBnd]));                          def = -def;                      }                      md->bonds[iBnd].z = def;                      cur_bnd = &(md->bondTypes[def]);                  }             }  // end if (save_lt)             // calculate distance (only once)             if (need_r)             {                dx = md->xyz[id1].x - md->xyz[id2].x;                dy = md->xyz[id1].y - md->xyz[id2].y;                dz = md->xyz[id1].z - md->xyz[id2].z;                delta_periodic(dx, dy, dz, md);                r2 = dx * dx + dy * dy + dz * dz;                need_r = 0;             }             action = 0;   // 0 - just cultivate bond 1 - delete bond 2 - transform bond             if ((cur_bnd->mxEx) && (r2 > cur_bnd->r2max))             {                 mnmx = 1;                 if (cur_bnd->new_type[mnmx] == 0)  // delete bond                   action = 1;                else                   action = 2;   // modify bond             }             else if ((cur_bnd->mnEx) && (r2 < cur_bnd->r2min))             {                 mnmx = 0;                 action = 2;   // at minimum only bond modification possible             }             // end select action             // try to change atom types (if needed)             if (action)             {                 save_lt = 1;                 new_spec1 = cur_bnd->new_spec1[mnmx];                 new_spec2 = cur_bnd->new_spec2[mnmx];                 //the first atom                 old = atomicCAS(&(md->types[id1]), spec1, new_spec1);                 if (old != spec1)                 {                     spec1 = old;                     spec2 = md->types[id2];   // refresh type of the 2nd atom                     // return to begin of the while loop                 }                 else      // types[id1] have been changed                 {#ifdef USE_NEWANG   // save changes in atom type                     atomicCAS(&(md->oldTypes[id1]), -1, spec1);#endif                     old_spec2 = spec2;                     while ((old = atomicCAS(&(md->types[id2]), old_spec2, new_spec2)) != old_spec2)                     {                         //! the worst variant: this thread changes atom 1, other thread changes atom 2                         // imagine that we had A-old bond with the same behavior                         def = md->def_bonds[spec1][old];#ifdef DEBUG_MODE                         if (def == 0)                         {                             printf("UBEH[001]: in apply_bonds, change atom types. There are no bond types between Species[%d] and [%d]\n", spec1, old);                             break;                         }#endif                         if (def < 0)  // spec1 -> new_spec2 spec2 -> newSpec1                         {                             cur_bnd = &(md->bondTypes[-def]);                             new_spec2 = cur_bnd->new_spec1[mnmx];                         }                         else // direct order                         {                             cur_bnd = &(md->bondTypes[def]);                             new_spec2 = cur_bnd->new_spec2[mnmx];                         }                         old_spec2 = old;#ifdef DEBUG_MODE                         cnt++;                         if (cnt > 10)                         {                             printf("UBEH[002]: too many atempst to change spec2 = %d\n", spec2);                             break;                         }#endif                     }#ifdef USE_NEWANG   // save changes in atom type                     atomicCAS(&(md->oldTypes[id2]), -1, spec2);#endif                     loop = 0;                 }                 //end change types             } // end if (action)             else               loop = 0;    // action == 0, out of cycle          }  // end while(loop)          if (action == 2)          {              new_bond_type = cur_bnd->new_type[mnmx];              if (new_bond_type < 0)              {                  md->bonds[iBnd].x = id2;                  md->bonds[iBnd].y = id1;                  new_bond_type = -new_bond_type;              }              md->bonds[iBnd].z = new_bond_type;              cur_bnd = &(md->bondTypes[new_bond_type]);          }          // perform calculations and save mean bond length          if (action != 1)  // not delete          {              r = sqrt(r2);              f = cur_bnd->force_eng(r2, r, eng, cur_bnd);              atomicAdd(&(md->frs[id1].x), f * dx);              atomicAdd(&(md->frs[id2].x), -f * dx);              atomicAdd(&(md->frs[id1].y), f * dy);              atomicAdd(&(md->frs[id2].y), -f * dy);              atomicAdd(&(md->frs[id1].z), f * dz);              atomicAdd(&(md->frs[id2].z), -f * dz);                            atomicAdd(&(cur_bnd->rSumm), r);              atomicAdd(&(cur_bnd->rCount), 1);          }          else      //delete bond          {// decrease the number of bonds for atoms              atomicSub(&(md->nbonds[id1]), 1);              atomicSub(&(md->nbonds[id2]), 1);              md->bonds[iBnd].z = 0;              // change parents              exclude_parents(id1, id2, md);          }          if (save_lt)          {              keep_bndlifetime(iStep, &(md->bonds[iBnd]), old_bnd);              if (action != 1) // not delete                atomicAdd(&(cur_bnd->count), 1);              atomicSub(&(old_bnd->count), 1);          }      } // end main loop      // split energy to shared and then to global memory      atomicAdd(&shEng, eng);      __syncthreads();      if (threadIdx.x == 0)          atomicAdd(&(md->engBond), shEng);}

В коде я использовал директивы препроцессора, чтобы включить проверки на те ситуации, которые могут возникнуть по недосмотру пользователя. Для ускорения быстродействия их можно отключить. Функция реализует указанную выше схему, но завернута ещё в один цикл, который пробегает по диапазону связей, за который отвечает данный поток. Здесь и далее идентификатор типа связи может быть отрицательный, это означает, что порядок атомов в связи нужно поменять местами (например, связь O-H и H-O это одна и та же связь, но в алгоритме порядок важен, чтобы это указать я использую индексы с противоположным знаком), делает это функция invert_bond, слишком тривиальная, чтобы её описывать. Функция delta_periodic применяет периодические граничные условия к разнице координат. Если нам нужно менять не только связи, но и валентные углы (директива USE_NEWANG), то необходимо помечать атомы, у которых мы поменяли тип (об этом далее). Для исключения повторного связывания одних и тех же атомов связью, массив parents хранит индекс одного из атомов связанных с данным (эта подстраховка работает не во всех случаях, но для моих этого достаточно). Если мы рвем какую-то связь, то нужно убрать из массива parents соответствующие индексы атомов, это делает функция exclude_parents:
__device__ void exclude_parents(int id1, int id2, cudaMD* md)// exclude id1 and id2 from parents of each other (if they are)//  and seek other parents if able{    // flags to clear parent    int clear_1 = 0;        int clear_2 = 0;    int i, flag;        if (md->parents[id1] == id2)         clear_1 = 1;    if (md->parents[id2] == id1)        clear_2 = 1;    i = 0;    while ((i < md->nBond) && (clear_1 || clear_2))    {        if (md->bonds[i].z != 0)        {            flag = 0;            if (clear_1)            {                if (md->bonds[i].x == id1)                {                    md->parents[id1] = md->bonds[i].y;                    flag = 1;                }                else if (md->bonds[i].y == id1)                {                    md->parents[id1] = md->bonds[i].y;                    flag = 1;                }                if (flag)                {                    clear_1 = 0;                    i++;                    continue;                }            }            if (clear_2)            {                if (md->bonds[i].x == id2)                {                    md->parents[id2] = md->bonds[i].y;                    flag = 1;                }                else if (md->bonds[i].y == id2)                {                    md->parents[id2] = md->bonds[i].y;                    flag = 1;                }                if (flag)                {                    clear_2 = 0;                    i++;                    continue;                }            }        }        i++;    }    // be on the safe side    if (clear_1)            md->parents[id1] = -1;    if (clear_2)        md->parents[id2] = -1;}

Функция, к сожалению, пробегает по всему массиву связей. Обрабатывать и удалять связи мы научились, теперь надо научится их создавать. Следующая функция помечает атомы, пригодные для образования химической связи:
__device__ void try_to_bind(float r2, int id1, int id2, int spec1, int spec2, cudaMD *md){    int r2Int;      //  (int)r2 * const    // save parents to exclude re-linking    if (md->parents[id1] == id2)        return;    if (md->parents[id2] == id1)        return;    if (md->bindBonds[spec1][spec2] != 0)    {        if (r2 < md->bindR2[spec1][spec2])        {            r2Int = (int)(r2 * 100);            if (atomicMin(&(md->r2Min[id1]), r2Int) > r2Int)    // replace was sucessfull            {                md->neighToBind[id1] = id2 + 1;     // as 0 is reserved for no neighbour                md->canBind[id1] = 1;            }            // similar for the second atom            if (atomicMin(&(md->r2Min[id2]), r2Int) > r2Int)    // replace was sucessfull            {                md->neighToBind[id2] = id1 + 1;     // as 0 is reserved for no bind                md->canBind[id2] = 1;            }        }    }}

В матрице bindBonds хранится информация о том, могут ли данные типы атомов образовывать связь и если могут то какую. В матрице bindR2 хранится максимальное расстояние между атомами необходимое для связывания. Если все условия благоприятные, то проверяем, нет ли у атомов других соседей пригодных для связывания, но поближе. Информация об ближайшем расстоянии до соседа хранится в массиве r2Min (для удобства массив имеет тип int и значения приводятся к нему с умножением на константу, 100). Если обнаруженный сосед наиближайший, то запоминаем его индекс в массиве neighToBind и ставим флаг canBind. Есть правда опасность, что пока мы перешли к обновлению индекса, другой поток перезаписал значение минимума, но это не так критично. Данную функцию целесообразно вызывать в функциях, выполняющих обход по парам атомов, например cell_list или all_pair, описанных в первой части. Само связывание:
__global__ void create_bonds(int iStep, int atPerBlock, int atPerThread, cudaMD* md)// connect atoms which are selected to form bonds{    int id1, id2, nei;    // neighbour index    int btype, bind;    // bond type index and bond index    cudaBond* bnd;    int spec1, spec2;   // species indexes        int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    int iat;    for (iat = id0; iat < N; iat++)    {        nei = md->neighToBind[iat];        if (nei)    // neighbour exists        {            nei--;  // (nei = spec_index + 1)            if (iat < nei)            {                id1 = iat;                id2 = nei;            }            else            {                id1 = nei;                id2 = iat;            }                        // try to lock the first atom            if (atomicCAS(&(md->canBind[id1]), 1, 0) == 0)  // the atom is already used                continue;            // try to lock the second atom            if (atomicCAS(&(md->canBind[id2]), 1, 0) == 0)  // the atom is already used            {                // unlock the first one back                atomicExch(&(md->canBind[id1]), 1);                continue;            }            // create bond iat-nei            bind = atomicAdd(&(md->nBond), 1);#ifdef DEBUG_MODE            if (bind >= md->mxBond)            {                printf("UBEH[003]: Exceed maximal number of bonds, %d\n", md->mxBond);            }#endif            spec1 = md->types[id1];            spec2 = md->types[id2];#ifdef USE_NEWANG   // save that we have changed atom type            atomicCAS(&(md->oldTypes[id1]), -1, spec1);            atomicCAS(&(md->oldTypes[id2]), -1, spec2);#endif            btype = md->bindBonds[spec1][spec2];                        if (btype < 0)            {                // invert atoms order                md->bonds[bind].x = id2;                md->bonds[bind].y = id1;                md->bonds[bind].z = -btype;                bnd = &(md->bondTypes[-btype]);                // change atom types according the formed bond                md->types[id1] = bnd->spec2;                md->types[id2] = bnd->spec1;            }            else             {                md->bonds[bind].x = id1;                md->bonds[bind].y = id2;                md->bonds[bind].z = btype;                bnd = &(md->bondTypes[btype]);                // change atom types according the formed bond                md->types[id1] = bnd->spec1;                md->types[id2] = bnd->spec2;            }                        atomicAdd((&bnd->count), 1);            md->bonds[bind].w = iStep;  // keep time of the bond creation for lifetime calculation            atomicAdd(&(md->nbonds[id1]), 1);            atomicAdd(&(md->nbonds[id2]), 1);            // replace parents if none:            atomicCAS(&(md->parents[id1]), -1, id2);            atomicCAS(&(md->parents[id2]), -1, id1);        }    }    // end loop by atoms}

Функция блокирует один атом, затем второй и, если удалось заблокировать оба, создает связь между ними. В самом начале функции идёт упорядочивание индексов атомов по порядку, чтобы исключить ситуацию, когда один поток блокирует первый атом в паре, а другой второй атом в той же паре, оба потока успешно проходят первую проверку и заваливаются на второй, в результате ни один из них не создает связь. Ну и, наконец, нам нужно удалять те связи, которые мы пометили на удаление в функции apply_bonds:
__global__ void clear_bonds(cudaMD* md)// clear bonds with .z == 0{    int i = 0;    int j = md->nBond - 1;    while (i < j)    {        while ((md->bonds[j].z == 0) && (j > i))            j--;        while ((md->bonds[i].z != 0) && (i < j))            i++;        if (i < j)        {            md->bonds[i] = md->bonds[j];            md->bonds[j].z = 0;            i++;            j--;        }    }    if ((i == j) && (md->bonds[i].z == 0))        md->nBond = j;    else        md->nBond = j + 1;}

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

Валентные углы


С валентными углами я введу несколько допущений, чтобы облегчить алгоритмы и функции, в результате они будут даже проще, чем для валентных связей. Во-первых, параметры валентных углов должны зависеть от всех трех атомов, но здесь мы будем считать, что тип валентного угла определяет исключительно атом в его вершине. Алгоритм образования/удаления углов предлагаю простейший: всегда, когда изменяли тип атома, запоминаем этот факт в соответствующем массиве oldTypes[]. Размер массива равен числу атомов, изначально он заполнен -1. Если какая-то функция меняет тип атома, она заменяет -1 на индекс первоначального типа. У всех атомов, которые поменяли тип, удаляем все валентные углы и пробегаем по всем связям этого атома, чтобы добавить соответствующие углы:
__global__ void refresh_angles(int iStep, int atPerBlock, int atPerThread, cudaMD *md)// delete old angles and create new ones for atoms which change their type{int i, j, n, t, ang;int nei[8];// bonded neighbors of given atomint cnt;int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;int N = min(id0 + atPerThread, md->nAt);int iat;for (iat = id0; iat < N; iat++)if (md->oldTypes[iat] != -1){i = 0;n = md->nangles[iat];while (n && (i < md->nAngle)){if (md->angles[i].w)if (md->angles[i].x == iat){n--;md->angles[i].w = 0;}i++;}// create new anglest = md->specAngles[md->types[iat]];// get type of angle, which formed by current atom typeif (t && (md->nbonds[iat] > 1))// atom type supports angle creating and number of bonds is enough{// search of neighbors by bondsi = 0; cnt = 0;n = md->nbonds[iat];while (n && (i < md->nBond)){if (md->bonds[i].z)// if bond isn't deleted{if (md->bonds[i].x == iat){nei[cnt] = md->bonds[i].y;cnt++;n--;}else if (md->bonds[i].y == iat){nei[cnt] = md->bonds[i].x;cnt++;n--;}}i++;}// add new angles based on found neighbors:for (i = 0; i < cnt-1; i++)for (j = i + 1; j < cnt; j++){ang = atomicAdd(&(md->nAngle), 1);md->angles[ang].x = iat;md->angles[ang].y = nei[i];md->angles[ang].z = nei[j];md->angles[ang].w = t;}n = (cnt * (cnt - 1)) / 2;}md->nangles[iat] = n;// reset flagmd->oldTypes[iat] = -1;}}

Массив specAngles содержит идентификаторы валентных углов, соответствующих данному типу атома. Следующая функция вызывает вычисление энергии и сил для всех углов:
__global__ void apply_angles(int iStep, int angPerBlock, int angPerThread, cudaMD* md)// apply valence angle potentials{cudaAngle* ang;// energies of angle potentialfloat eng;__shared__ float shEng;if (threadIdx.x == 0)shEng = 0.0f;__syncthreads();int id0 = blockIdx.x * angPerBlock + threadIdx.x * angPerThread;int N = min(id0 + angPerThread, md->nAngle);int i;for (i = id0; i < N; i++)if (md->angles[i].w){ang = &(md->angleTypes[md->angles[i].w]);ang->force_eng(&(md->angles[i]), ang, md, eng);}// split energy to shared and then to global memoryatomicAdd(&shEng, eng);__syncthreads();if (threadIdx.x == 0)atomicAdd(&(md->engAngl), shEng);}

Ну и для примера потенциала таких углов, даю функцию гармонического косинуса, на которую может указывать поле force_eng структуры cudaAngle:
__device__ void angle_hcos(int4* angle, cudaAngle* type, cudaMD* md, float& eng)// harmonic cosine valent angle potential:// U = k / 2 * (cos(th)-cos(th0))^{float k = type->p0;float cos0 = type->p1;// indexes of central atom and ligands:int c = angle->x;int l1 = angle->y;int l2 = angle->z;// vector ijfloat xij = md->xyz[l1].x - md->xyz[c].x;float yij = md->xyz[l1].y - md->xyz[c].y;float zij = md->xyz[l1].z - md->xyz[c].z;delta_periodic(xij, yij, zij, md);float r2ij = xij * xij + yij * yij + zij * zij;float rij = sqrt(r2ij);// vector ikfloat xik = md->xyz[l2].x - md->xyz[c].x;float yik = md->xyz[l2].y - md->xyz[c].y;float zik = md->xyz[l2].z - md->xyz[c].z;delta_periodic(xik, yik, zik, md);float r2ik = xik * xik + yik * yik + zik * zik;float rik = sqrt(r2ik);float cos_th = (xij * xik + yij * yik + zij * zik) / rij / rik;float dCos = cos_th - cos0; // delta cosinusfloat c1 = -k * dCos;float c2 = 1.0 / rij / rik;atomicAdd(&(md->frs[c].x), -c1 * (xik * c2 + xij * c2 - cos_th * (xij / r2ij + xik / r2ik)));atomicAdd(&(md->frs[c].y), -c1 * (yik * c2 + yij * c2 - cos_th * (yij / r2ij + yik / r2ik)));atomicAdd(&(md->frs[c].z), -c1 * (zik * c2 + zij * c2 - cos_th * (zij / r2ij + zik / r2ik)));atomicAdd(&(md->frs[l1].x), c1 * (xik * c2 - cos_th * xij / r2ij));atomicAdd(&(md->frs[l1].y), c1 * (yik * c2 - cos_th * yij / r2ij));atomicAdd(&(md->frs[l1].z), c1 * (zik * c2 - cos_th * zij / r2ij));atomicAdd(&(md->frs[l2].x), c1 * (xij * c2 - cos_th * xik / r2ik));atomicAdd(&(md->frs[l2].y), c1 * (yij * c2 - cos_th * yik / r2ik));atomicAdd(&(md->frs[l2].z), c1 * (zij * c2 - cos_th * zik / r2ik));eng += 0.5 * k * dCos * dCos;}

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

Примеры


Не претендуя на точность, я попробовал изобразить сборку молекул воды из отдельных ионов. Парные потенциалы задал произвольно в форме потенциала Букингема, а затем добавил возможность создавать связи в виде гармонического потенциала, с равновесным расстоянием равным длине связи H-O в воде, 0.96 . Кроме того, при связывании второго протона с кислородом добавлялся валентный угол с вершиной в кислороде. После 100000 шагов из случайно разбросанных ионов появились первые молекулы. На рисунке показаны начальная (слева) и конечная (справа) конфигурации:

Можно поставить вот такой эксперимент: пусть вначале все атомы одинаковые, однако оказавшись рядом с друг другом они образуют связь. Позволим связанным атомам образовывать ещё одну связь либо со свободным атомом, либо с ещё одной такой же связанной молекулой. В результате у нас получится вот какая-то такая самоорганизация, где атомы выстраиваются в цепочки:


Заключительные комментарии


1. Здесь для связывания мы использовали лишь один критерий расстояние, хотя в принципе, могут быть и другие, например энергия системы. В реальности при образовании химической связи, как правило, энергия выделяется в виде теплоты. Здесь это никак не учтено, но можно попробовать это реализовать, например изменять скорости частиц.
2. Взаимодействия между частицами через потенциал химической связи не отменяют того, что частицы могут ещё взаимодействовать через межмолекулярные парные потенциалы и Кулоновское взаимодействие. Можно было бы, конечно, для связанных атомов не вычислять межмолекулярные взаимодействия, но это, в общем случае, требует долгих проверок. Поэтому проще задать потенциал химической связи таким образом, чтобы его сумма с другими потенциалами давала нужную функцию.
3. Параллельная реализация связывания частиц не только дает выигрыш в скорости, но даже выглядит более реалистично, поскольку частицы соревнуются между собой.

Ну и на хабре есть несколько проектов, очень близких к данному:
Подробнее..

Что такое температура и как её учитывать в молекулярном моделировании? Реализация на CUDA

06.01.2021 20:05:30 | Автор: admin

Методы численного моделирования молекулярных систем, такие как молекулярная динамика, рассматривают эти системы как механические (что-то вроде набора шариков на пружинках). Однако, в отличие от механических систем, для молекулярных существует понятие температура. Вещество не может существовать без температуры, а температура без вещества (на счет последней части утверждения есть и другое мнение). Из опыта мы знаем, что очень многие свойства вещества кардинально зависят от температуры, и, естественно, что её надо как-то учитывать. Для поддержания температуры в молекулярной динамике используются специальные алгоритмы-термостаты. Наиболее известные среди них это термостаты Андерсена, Берендсена и Нозе-Гувера. Все они основаны на молекулярно-кинетической теории газов, где температура есть просто величина пропорциональная среднекинетической энергии молекул. Соответственно, работа данных термостатов осуществляется путём умножения скоростей частиц на некоторую величину.

Что же не так?

  1. Масштабирование скоростей изменяет их модуль, но не направление. А ведь тепловое движение, согласно большинству утверждений, должно быть хаотическим. Оно определяет броуновское движение. Термостатирование в молекулярной динамике же не привносит никакого хаоса, частицы движутся в том же направлении, как и без учета температуры, только быстрее/медленнее. Метод остается полностью детерминированным, начальные координаты и скорости частиц полностью определяют конечное состояние.

  2. Скорости определяются через интегрирование уравнений движения, из ускорений а ускорения из сил. Силы задаются потенциалами, которые позиционируются как имеющие смысл. Получается мы сначала выбрали потенциал, через него определили как частицы взаимодействуют, начали вычислять ускорения и скорости, а потом взяли эти скорости и домножили на какую-то величину. Какой тогда смысл в потенциалах?

  3. Связь температуры и кинетической энергии выводится из основных положений МКТ, которые верны только для идеального газа. Наиболее близкими к идеальному газу реальные системы это благородные газы при высоких температурах и низких давлениях. Проиллюстрировать неуниверсальность МКТ можно простым примером: при одной температуре средняя кинетическая энергия атомов в твердом теле и в газе должна быть одинаковой, что не кажется правдоподобным.

Модель

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

Соберем известные факты о температуре:

  • Температура не может быть ниже абсолютного нуля, т.е. 0 К.

  • Температура понятие интенсивное, оно не зависит от размеров системы, количества вещества и т.д.

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

  • Скорости химических реакций и коэффициенты диффузии с температурой растут экспоненциально.

  • При повышении температурах сложные вещества постепенно разлагаются на более простые.

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

  • С увеличением температуры тела обычно расширяются (газы все, твердые тела почти все, но есть исключения, см. отрицательный КТР).

  • Теплоемкость с температурой как правило растет (для некоторых жидкостей, такие как расплавленные оксиды она постоянная или чуть падает).

  • С температурой связан интересный эффект, адиабатическое расширение в пустоту (эффект Джоуля-Томсона).

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

  • С увеличением температуры растет сжимаемость, т.е. при высоких температурах чтобы сжать тело в то же число раз, нужно приложить меньшее давление. По крайней мере, это верно для расплавов.

  • Электропроводность металлов падает с температурой.

  • Спектр излучения абсолютно черного тела определяется исключительно его температурой.

Как мы видим, количественных характеристик температуры не так уж и много. Привязываться к каким-то ключевым точкам типа температуры кипения воды не имеет смысла, поскольку: они очень частные; зависят и от других условий, в первую очередь, давления. Тем не менее, есть характеристика, кажущаяся довольно перспективной. Это спектр излучения абсолютно черного тела. Оказывается, если тело ничего не отражает, а все излучение от него является собственным, то спектр этого излучения не зависит от природы объекта, а определяется исключительно его температурой, согласно формуле Планка. Тот факт, что все тела излучают и характеристики этого излучения зависят от температуры используется, например, в оптической пирометрии. Таким образом, температуру можно определить через спектр излучения. В условии равновесия с термостатом, т.е. средой, имеющей заданную температуру, излучение, испускаемое объектом, должно сравняться с излучением термостата. Тогда можно реализовать термостатирование следующим образом: атомы поглощают и испускают фотоны в случайном направлении, при этом их скорости изменяются согласно закону сохранения импульса. Интуитивно кажется, что если мы будем прибавлять/убавлять поправку к скорости атома в случайном направлении, то атом сначала наберет какую-то скорость, а затем эта скорость выйдет на насыщение, поскольку все выше вероятность того, что новая скорость будет направлена против одной из предыдущих добавок. Однако легко показать, что это не так, начнем с закона сохранения импульса:

m\vec{v_1}=m\vec{v_0}+\vec{p_\gamma},\qquad\qquad\qquad (1)

где m и v масса и скорость атома, индексы 1 и 0 означают состояния до и после поглощения фотона, p импульс фотона. Перейдем к скоростям:

\vec{v_1}=\vec{v_0}+\frac\epsilon{mc}\vec{r},\quad\quad\quad(2)

где энергия фотона, c скорость света в вакууме, r единичный вектор случайного направления. Из этого, квадрат скорости равен:

v_1^2 = (\vec{v_0} + \frac{mc}\vec{r})^2 = v_0^2 + 2 \frac{mc}\vec{r}\vec{v_0} + \frac{^2}{m^2c^2}, \quad\quad\quad(3)

или, используя определение скалярного произведения векторов:

v_1^2 = v_0^2 + 2 \frac{mc}v_0 \cos \phi + \frac{^2}{m^2c^2},\quad\quad(4)

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

\langle\Delta v^2\rangle = \frac{^2}{m^2c^2}\quad\quad(5)

положительное. Т.е. действуя таким образом, мы будем накачивать систему кинетической энергией, независимо от того будем излучать или поглощать фотоны, поскольку знак будет стоять перед косинусом, который в среднем нулевой. Чтобы скомпенсировать неограниченный рост скоростей в системе нужно, чтобы среднее изменение квадрата скорости атома при излучении фотона было равно среднему значению изменения квадрата скорости при поглощении фотона (5) с обратным знаком:

\langle\Delta v^2\rangle = -\frac{^2}{m^2c^2}\quad\quad(6)

для этого косинус между направлением движения атома и излучаемым фотоном, в среднем, должен равняться:

\langle\cos\phi\rangle = -\frac{}{2mcv_0}\quad\quad(7)

Это выражение устанавливает связь между скоростями атомов и энергией фотонов. Можно задать косинус равномерной случайной величиной, так чтобы указанное выше значение было серединой интервала, а -1 его нижней границей.

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

\frac12mv_1^2 = \frac12mv_0^2 + ,\quad\quad(8)

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

U_1 + K_1 = U_0 + K_0 + ,\quad\quad(9)

где K кинетическая энергия. Остался вопрос, фотоны какой энергии излучать атому? Давайте немного упростим модель, пусть тепловое излучение описывается не формулой Планка, а будет монохроматическим, с длиной волны / частотой / энергией однозначно определяемыми температурой. Для этого подходит, например, закон Вина, определяющий длину волны, дающую наибольший вклад в спектр излучения:

[м] = 0.002898 / T,\quad\quad(10)

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

 = T,\quad\quad(11)

где некоторый постоянный коэффициент. Т.е. температура T, это такое состояние вещества, когда излучаются фотоны энергией T.

Кроме того, у нас есть ещё одна численная характеристика, касающаяся температуры. Мы совершенно точно знаем, какое количество теплоты (энергии/тепловой энергии) нужно сообщить телу, чтобы нагреть его до заданной температуры от нуля кельвинов. В условиях постоянного давления эта величина называется приращением энтальпии / инкрементом энтальпии / энтальпией (не путать с энтальпией образования), её значение зависит не только от температуры, но и от природы вещества. Если приводить образную аналогию, то вещество можно представить как сосуды с некой жидкостью, где её уровень это температура, а объём энтальпия. Производная энтальпии по температуре называется изобарной теплоемкостью, cp. Для идеального газа (и, с некоторым приближением, для благородных газов) cp постоянна, тогда:

H = c_p T.\quad\quad(12)

Объединяя (11) и (12) можно исключить температуру и получить энергию фотона как функцию от текущей тепловой энергии, запасенной в веществе:

 = \frac{c_p} H.\quad\quad(13)

Мы получили, что энергия излучаемых фотонов определяется энтальпией вещества. Продолжая нашу аналогию с сосудами, можно представить, что они теряют жидкость пропорционально высоте столба. Напоминаю, что выражение (13) справедливо только для систем с постоянной теплоемкостью, для реальных систем нужно находить температуру как функцию от энтальпии и подставлять её в (11). При переходе к числам меня ждал сюрприз, выразив из (10) значение и подставив его в (13), а cp взяв равным 5/2k (изобарная теплоемкость одноатомного идеального газа) получилось, что коэффициент перед энтальпией больше 1, а это означает, что система должна излучать больше энергии, чем имеется, что, конечно, не физично. Немного покопавшись, я нашел, что при переходе от длины волны к частоте (и затем к энергии) нужно использовать несколько отличное выражение, чем (10). В конечном итоге, с использованием данных по теплоемкости гелия, у меня получилось, /cp 0.9, что я и взял для своих численных экспериментов. В качестве энтальпии я взял величину U, описанную выше, формула (9).

Немного кода

Основная часть кода базируется на изложенном ранее. Однако вместо термостата Нозе-Гувера используется вышеописанная схема. Функция термостатирования:

__global__ void tstat_radi(int atPerBlock, int atPerThread, cudaMD* md){    __shared__ int indEng;    __shared__ float engTemp;    if (threadIdx.x == 0)    {        indEng = atomicAdd(&(md->curEng), 1);   // get current index in photon energies array        engTemp = 0.f;    }    __syncthreads();    uint4 randVar = md->ui4rnd;    float pe;        // photon energy    int i, e0, ei;// indexes    // atomic range    int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    e0 = indEng * atPerBlock + threadIdx.x * atPerThread - id0;    for (i = id0; i < N; i++)    {        ei = e0 + i;        if (ei >= md->nAt)            ei -= md->nAt;        pe = md->engPhotons[ei];        adsorb_rand_photon(&(md->vls[i]), &(md->engs[i]), md->masses[i], pe, randVar, md);        radiate_photon(&(md->vls[i]), &(md->engs[i]), md->masses[i], randVar, md);    }    atomicAdd(&engTemp, teng);    if (threadIdx.x == 0)    {        rnd_xor128(md->ui4rnd);        atomicAdd(&(md->engTemp), engTemp);    }}

Функция принимает на вход число атомов в потоке и блоке, а также указатель на структуру, где хранятся все МД-данные. Для каждого атома в обрабатываемом диапазоне находится энергия случайного фотона, вызываются функции его поглощения и функция излучения. Энергия поглощаемого фотона берется из предварительно сгенерированного массива энергий, которое соответствует Планковскому распределению (эта часть вычисляется на хосте и загружается в массив). Я здесь использовал ту фичу CUDA, что порядок выполнения блоков не определён. Первый запустившийся блок берет первый кусок массива энергий фотонов и т.д., таким образом получается что-то вроде рандомизации. Таким же образом я хотел сделать и направления фотонов, но рандомизация выходит слабенькой, поскольку порядок выполнения блоков хоть и не определен, но все-таки даёт очень небольшой разбег, для направлений это критично, даёт заметные артефакты.

Функция поглощения фотона принимает на вход указатель на скорость атома, указатель на тепловую энергию атома / энтальпию, массу атома, энергию фотона, переменную для генератора случайных чисел и указатель на главную структуру:

__device__ void adsorb_rand_photon(float3 *vel, float *int_eng, float mass, float eng, uint4 &rand_var, cudaMD *md){    float u0 = *int_eng;    float v02 = float3_sqr(*vel);   // square of initial velocity    float3 rand_vect = rand_uvect(rand_var, md);    // momentum conservation:    float ermc = eng * revLight / mass;    vel->x += ermc * rand_vect.x;     vel->y += ermc * rand_vect.y;    vel->z += ermc * rand_vect.z;    float v12 = float3_sqr(*vel);    // energy conservation: old kinetic energy + photon energy = new kinetic energy + 'internal' energy    *int_eng += eng + 0.5f * mass * (v02 - v12);}

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

__device__ void radiate_photon(float3* vel, float* int_eng, float mass, uint4& rand_var, cudaMD* md){    float u0 = *int_eng;    float v02 = float3_sqr(*vel);   // square of initial velocity    float v0 = sqrt(v02);           // module of initial velocity    float ph_eng = 0.9f * u0;    float ermc = ph_eng * revLight / mass;    // random radiation    float3 rand_vect;    // random cosine between   (1-2a/v0) .. -1 with mean = -a/v0, where a = e/mc    float ermcv0 = ermc / v0;       // ph_eng/(m*c*v0)    float cos_phi;    if (ermcv0 >= 1.f)        cos_phi = -1.f;    else    {        int rnd = rnd_xor128(rand_var) % 2048;        float cos_phi = (float)rnd / 1024.f * (1.f - ermcv0); // 1024, because factor of 2          cos_phi -= 1.f;        rnd = rnd_xor128(rand_var) % 2048;        float theta = (float)rnd / 1024.f * numPi;  // 0 .. 2PI        rand_vect = get_angled_vector(*vel, cos_phi, theta);    }    // momentum conservation:    vel->x += ermc * rand_vect.x;       vel->y += ermc * rand_vect.y;    vel->z += ermc * rand_vect.z;    float v12 = float3_sqr(*vel);    // energy conservation:    *int_eng -= (ph_eng + 0.5f * mass * (v12 - v02));}

Тут учитывается, что выражение (7) может дать и косинус меньше -1, что нужно пофиксить. Поскольку в 3-мерном пространстве векторов, находящихся под заданным углом к данному, целый конус, функция get_angled_vector вычисляет один из них, по выбранному углу вращения:

__device__ float3 get_angled_vector(float3 invec, float cos_phi, float theta){    float3 v1, v2, v3;    float leng1 = float3_length(invec);    // only for non-zero input vectors!    v1 = make_float3(invec.x / leng1, invec.y / leng1, invec.z / leng1);    // find v2 perpendicular to v1    if (v1.x != 0.f)    {        v2.y = 1.f; v2.z = 1.f; // any coordinates        v2.x = -(v1.y * v2.y + v1.z * v2.z) / v1.x;    }    else        if (v1.y != 0.f)        {            v2.x = 1.f; v2.z = 1.f; // any coordinates            v2.y = -(v1.z * v2.z) / v1.y;     // a1 = 0 !        }        else // a1=0, b1=0, c1 <> 0:        {            v2.x = 1.f; v2.y = 0.f; v2.z = 0.f;        }    // v3 is perpendicular to both v1 and v2    v3.x = v1.y * v2.z - v1.z * v2.y;    v3.y = -v1.x * v2.z + v1.z * v2.x;    v3.z = v1.x * v2.y - v1.y * v2.x;    // convert v2 and v3 into unit vectors    float leng2 = float3_length(v2);    float leng3 = float3_length(v3);    v2.x /= leng2;    v2.y /= leng2;    v2.z /= leng2;    v3.x /= leng3;    v3.y /= leng3;    v3.z /= leng3;    float sinPhi, sinTh, cosTh;    sinPhi = sqrt(1 - cos_phi * cos_phi);    sincos(theta, &sinTh, &cosTh);    // parameteric equation of circle    float x = v1.x * cos_phi + sinPhi * (cosTh * v2.x + sinTh * v3.x);    float y = v1.y * cos_phi + sinPhi * (cosTh * v2.y + sinTh * v3.y);    float z = v1.z * cos_phi + sinPhi * (cosTh * v2.z + sinTh * v3.z);    float3 res = make_float3(x, y, z);    return res;} 

Генератор случайных чисел на CUDA не тривиальная задача, спасибо пользователю @denglide за обзор генераторов и пользователю @maratyszcza за xorShift генератор, который и был использован в данной работе.

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

Результаты

Ну а теперь к численным экспериментам. В качестве модельного объекта я взял 40000 атомов аргона, взаимодействующих согласно парному потенциалу из этой работы. Размеры бокса выбраны таким образом, чтобы плотность системы соответствовали плотности аргона при стандартных условиях (T = 298 К и P = 1 атм). Все приведенные результаты даны для 1 млн шагов интегрирования с шагом интегрирования 1 фс. Вычисления проводил на видеокарте GeForce RTX 2080 Ti.

Рис. 1. Зависимость кинетической энергии от времени эволюции системы

Ещё раз напоминаю, что целевую температуру я задаю через набор энергий фотонов, излучаемых термостатом. На данном этапе я следил за изменением кинетической энергии системы (рис. 1) при различных температурах термостата и разных вариантов начальных скоростей атомов. В ряде случаев начальные скорости были нулевыми (Ekin = 0), в других случаях были для всех атомов одинаковыми по модулю, но направленными хаотично (Ekin > 0). Независимо от величины начальных скоростей, кинетическая энергия выходила на примерно постоянное значение, зависящее от температуры термостата, т.е. такой термостат умеет разгонять или замедлять атомы, в зависимости от текущего состояния системы и целевой температуры, что в первую очередь от него и требуется.

Рис. 2. Распределение атомов по скоростям после 1млн шагов эволюции системы.

Второе, на что я обратил внимание, это распределение скоростей атомов, после 1 млн шагов, рис. 2. Несмотря на то, что начальные скорости были заданы одинаковыми по модулю, в результате работы термостата распределение атомов по скоростям приобрело вид, очень похожий на распределение Максвелла. Точно также как и в распределении Максвелла, меньшая температура соответствует более узкому пику. То, что такое распределение получилось само собой, при выполнении простых правил, я считаю определенным успехом.

Заключение

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

Плюсы предложенного термостата:

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

  2. Распределение скоростей частиц в системе приобретает вид распределения Максвелла.

  3. Вычислительная сложность порядка O(N), как у других термостатов.

  4. Температура определяет скорости частиц, но не сводится к ним.

  5. Есть возможность учитывать индивидуальные свойства разных веществ, используя их температурную зависимость энтальпии (вернее энтальпийную зависимость температуры).

Минусы и недоработки:

  1. Сумма скоростей системы не равна нулю, возможно проблема в генераторе случайных чисел.

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

  3. Формула (7) выведена исключительно из математических соображений и имеет слабое физическое обоснование.

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

Ну и напоследок, вот как выглядят 40000 атомов:

Подробнее..

Категории

Последние комментарии

  • Имя: Макс
    24.08.2022 | 11:28
    Я разраб в IT компании, работаю на арбитражную команду. Мы работаем с приламы и сайтами, при работе замечаются постоянные баны и лаги. Пацаны посоветовали сервис по анализу исходного кода,https://app Подробнее..
  • Имя: 9055410337
    20.08.2022 | 17:41
    поможем пишите в телеграм Подробнее..
  • Имя: sabbat
    17.08.2022 | 20:42
    Охренеть.. это просто шикарная статья, феноменально круто. Большое спасибо за разбор! Надеюсь как-нибудь с тобой связаться для обсуждений чего-либо) Подробнее..
  • Имя: Мария
    09.08.2022 | 14:44
    Добрый день. Если обладаете такой информацией, то подскажите, пожалуйста, где можно найти много-много материала по Yggdrasil и его уязвимостях для написания диплома? Благодарю. Подробнее..
© 2006-2024, personeltest.ru