Method Article
Атмосферные концентрации слабо связанных молекулярных скоплений можно вычислить из термохимических свойств низкоэнергетических структур, найденных с помощью многоступенчатой методологии конфигурации выборки с использованием генетического алгоритма и полуэмпирической и ab initio квантовой химии.
Вычислительное исследование образования и роста атмосферных аэрозолей требует точной поверхности свободной энергии Гиббса, которая может быть получена из газовой фазы электронной структуры и вибрационных вычислений частоты. Эти количества действительны для тех атмосферных скоплений, геометрия которых соответствует минимуму по их потенциальным энергетическим поверхностям. Свободная энергия Гиббса минимальной энергетической структуры может быть использована для прогнозирования атмосферных концентраций скопления в различных условиях, таких как температура и давление. Мы представляем вычислительно недорогой процедуры, построенные на генетическом алгоритме на основе конфигурации выборки следуют ряд все более точные расчеты скрининга. Процедура начинается с генерации и развития геометрии большого набора конфигураций с использованием полуэмпирических моделей, а затем уточняет полученные уникальные структуры на ряде уровней теории ab initio высокого уровня. Наконец, термодинамические коррекции вычисляются для полученного набора минимальных энергетических структур и используются для вычисления свободных энергий образования, равновесных констант и атмосферных концентраций. Мы представляем применение этой процедуры для изучения гидратированных глицинных кластеров в условиях окружающей среды.
Наиболее неопределенным параметром в атмосферных исследованиях изменения климата является точная степень, в которой частицы облаков отражают поступающую солнечную радиацию. Аэрозоли, которые являются твердыми частицами, взвешенными в газе, образуют частицы облаков, называемые ядрами конденсации облаков (CCN), которые рассеивают поступающее излучение, тем самым предотвращая его поглощение и последующее нагревание атмосферы1. Детальное понимание этого эффекта чистого охлаждения требует понимания роста аэрозолей в CCNs, что, в свою очередь, требует понимания роста малых молекулярных скоплений в аэрозольные частицы. Недавняя работа показала, что образование аэрозолей инициируется молекулярными кластерами диаметром 3 нм или менее2; однако, этот режим размера трудно получить доступ с помощью экспериментальных методов3,4. Поэтому для преодоления этого экспериментального ограничения желательно принять подход к вычислительному моделированию.
Используя наш подход моделирования, описанный ниже, мы можем проанализировать рост любого гидратного кластера. Поскольку мы заинтересованы в роли воды в формировании крупных биологических молекул из небольших компонентов в добиотической среде, мы иллюстрируем наш подход глицином. Проблемы и инструменты, необходимые для решения этих вопросов исследования очень похожи на те, которые участвуют в изучении атмосферных аэрозолей и донуклеационных кластеров5,,6,,7,,8,9,,10,,11,12,13,14,15. Здесь мы исследуем гидратированные глицинные кластеры, начиная с изолированной молекулы глицина, за которыми следует серия пошаговых дополнений до пяти молекул воды. Конечная цель состоит в том, чтобы рассчитать равновесные концентрации скоплений Gly (H2O)n'0-5 в атмосфере при комнатной температуре на уровне моря и относительной влажности (RH) в 100%.
Небольшое число этих субнанометров молекулярных кластеров вырастают в метастабильное критическое скопление (1-3 нм в диаметре) либо путем добавления других молекул пара, либо свогуляющих на существующих кластерах. Эти критические кластеры имеют благоприятный профиль роста, ведущий к образованию гораздо больших (до 50-100 нм) ядер конденсации облаков (CCN), которые непосредственно влияют на эффективность осадков облаков, а также их способность отражать свет инцидента. Поэтому хорошее понимание термодинамики молекулярных скоплений и их равновесного распределения должно привести к более точным прогнозам воздействия аэрозолей на глобальный климат.
Описательная модель формирования аэрозолей требует точной термодинамики формирования молекулярного кластера. Вычисление точной термодинамики формирования молекулярного кластера требует идентификации наиболее стабильных конфигураций, что предполагает поиск глобальной и локальной минимальной на потенциальной энергетической поверхности кластера (PES)16. Этот процесс называется конфигурационной выборки и может быть достигнуто с помощью различных методов, в том числе на основе молекулярной динамики (MD)17,18,19,20, Монте-Карло (MC)21,22, и генетические алгоритмы (GA)23,24,25.
Различные протоколы были разработаны на протяжении многих лет для получения структуры и термодинамики атмосферных гидратов на высоком уровне теории. Эти протоколы отличались выбором (i) метода конфигурационной выборки, ii) характером метода низкого уровня, используемого в конфигурационной выборке, и iii) иерархией методов более высокого уровня, используемых для уточнения результатов в последующих шагах.
Конфигурационные методы отбора проб включали химическую интуицию26, случайную выборку27,,28,молекулярной динамики (MD)29,30, бассейн прыжков (BH)31, и генетический алгоритм (GA)24,25,32. Наиболее распространенными низкоуровневыми методами, используемыми с помощью этих методов отбора проб, являются силовые поля или полуэмпирические модели, такие как PM6, PM7 и SCC-DFTB. За ними часто следуют расчеты DFT с все более большими базовыми наборами и более надежными функциональными возможностями с более высоких ступеней лестницыИакова 33. В некоторых случаях за ними следуют методы волной более высокого уровня, такие как MP2, CCSD (T), а также экономически эффективные DLPNO-CCSD (T)34,35.
Kildgaard et al.36 разработали систематический метод, при котором молекулы воды добавляются в точках на сферах Фибоначчи37 вокруг небольших гидратированных или негидратированных кластеров для генерации кандидатов для больших кластеров. Нефизические и избыточные кандидаты удаляются на основе пороговых значений контакта и расстояния корня-квадрата между различными конформистами. Последующие оптимизации с использованием полуэмпирического метода PM6 и иерархии методов DFT и wavefunction используются для получения набора низкоэнергетических конформистов на высоком уровне теории.
Алгоритм искусственной пчелообразной колонии (ABC)38 представляет собой новый подход к выборке конфигурации, который недавно был реализован Чжан и др. для изучения молекулярных кластеров в программе под названием ABCluster39. Kubecka et al.40 использовали ABCluster для конфигурационной выборки с последующим низкоуровневым переоптимизацией с использованием герметичного полуэмпирического метода GFN-xTB41. Они дополнительно уточнили структуры и энергии с помощью методов DFT, а затем конечных энергий с использованием DLPNO-CCSD (T).
Независимо от метода, конфигурационная выборка начинается с случайно- или неслучайно генерируемого распределения точек на PES. Каждая точка соответствует конкретной геометрии рассматриваемого молекулярного кластера и генерируется методом выборки. Затем ближайший локальный минимум находится для каждой точки, следуя направлению «вниз» на PES. Таким образом, найденный набор минимы соответствует геометриям молекулярного скопления, которые стабильны, по крайней мере на некоторое время. Здесь форма PES и оценка энергии в каждой точке на поверхности будут чувствительны к физическому описанию системы, где более точное физическое описание приводит к более вычислительно-дорогому расчету энергии. Мы специально будем использовать метод GA, реализованный в программе OGOLEM25, который был успешно применен к различным глобальным проблемам оптимизации и конфигурации выборки42,,43,,44,,45, для создания первоначального набора точек отбора проб. PES будет описано PM7 модель46 реализованы в программе MOPAC201647. Эта комбинация используется потому, что она генерирует большее разнообразие точек по сравнению с методами MD и MC и находит локальные минимумы быстрее, чем более подробные описания PES.
Набор оптимизированных ГА локальных минимумов берется в качестве стартовой геометрии для серии скрининговых шагов, которые приводят к набору низколежащей минимальной энергии. Эта часть протокола начинается с оптимизации набора уникальных ОПтимизированных га структур с использованием плотной функциональной теории (DFT) с небольшим базовым набором. Этот набор оптимизаций, как правило, дает меньший набор уникальных локальных минимальных структур, которые моделируются более подробно по сравнению с ОПтимизированными ГА полуэмпирическими структурами. Затем на этом уменьшенном наборе структур выполняется еще один раунд оптимизации DFT с использованием более широкого базисного набора. Опять же, этот шаг, как правило, дают меньший набор уникальных структур, которые моделируются более подробно по сравнению с небольшой базой DFT шаг. Окончательный набор уникальных структур затем оптимизирован до более жесткой конвергенции и рассчитываются гармонические вибрационные частоты. После этого шага у нас есть все необходимое для вычисления равновесных концентраций скоплений в атмосфере. Общий подход обобщен диаграммонана на рисунке 1. Мы будем использовать PW9148 обобщенный градиент приближения (GGA) обмен-корреляции функциональных в Gaussian0949 реализации DFT вместе с двумя вариациями Pople50 базовый набор (6-31"G " для небольшой базовый шаг и 6-311 "G" для большого базисного шага). Данная комбинация валютно-корреляционного функционального и базисного набора была выбрана благодаря своему предыдущему успеху в вычислении точных свободных энергий образования Гиббса для атмосферных кластеров51,52.
Этот протокол предполагает, что пользователь имеет доступ к высокопроизводительным вычислительным (HPC) кластеру с портативной пакетной системой53 (PBS), MOPAC2016 (http://openmopac.net/MOPAC2016.html)47,OGOLEM (https://www.ogolem.org)25, Gaussian 09 (https://gaussian.com)49, и OpenBabel54 (http://openbabel.org/wiki/Main_Page программное обеспечение, установленное в соответствии с их конкретными инструкциями по установке. Каждый шаг в этом протоколе также использует набор внутренних оболочек и скриптов Python 2.7, которые должны быть сохранены в каталоге, который включен в $PATH экологической переменной пользователя. Все необходимые экологические модули и разрешения на выполнение всех вышеперечисленных программ также должны быть загружены в сеанс пользователя. Использование диска и памяти кодом GA (OGOLEM) и полуэмпирическими кодами (MOPAC) очень мало по современным стандартам компьютерных ресурсов. Общее использование памяти и диска для OGOLEM/MOPAC зависит от того, сколько потоков человек хочет использовать, и даже в этом случае использование ресурсов будет небольшим по сравнению с возможностями большинства систем HPC. Потребности в ресурсах методов КМ зависят от размера кластеров и уровня используемой теории. Преимущество использования этого протокола заключается в том, что можно варьировать уровень теории, чтобы иметь возможность вычислить окончательный набор низкоэнергетических структур, имея в виду, что обычно более быстрые расчеты приводят к большей неопределенности в точности результатов.
Для ясности локальный компьютер пользователя будет называться«локальным компьютером»,в то время как кластер HPC, к которому они имеют доступ, будет называться«удаленным кластером».
1. Поиск минимальной энергетической структуры изолированных глицина и воды
ПРИМЕЧАНИЕ: Цель здесь состоит в двух раз: i) получить минимальные энергетические структуры изолированных молекул воды и глицина для использования в генетическом алгоритме конфигурационных выборок, (ii) и вычислить термодинамические коррекции газовых фазовых энергий этих молекул для использования при расчете атмосферных концентраций.
2. Генетическо-алгоритм-на основе конфигурационной выборки Gly (H2O)n-1-5 кластеров
ПРИМЕЧАНИЕ: Цель здесь состоит в том, чтобы получить набор низкоэнергетических структур для Gly (H2O)n'1-5 на недорогом полуэмпирическом уровне теории, используя модель PM746, реализованную в MOPAC47. Крайне важно, чтобы рабочий каталог имеет точную организацию и структуру, как показано на рисунке 2. Это необходимо для обеспечения того, чтобы пользовательские скрипты оболочки и Python работали без сбоев.
3. Уточнение с использованием метода ЗМ с небольшим набором основы
ПРИМЕЧАНИЕ: Цель здесь заключается в уточнении конфигурационной выборки кластеров Gly (H2O)n-1-5 с использованием лучшего квантово-механического описания для получения меньшего, но более точного набора кластерных структур Gly (H2O)n-1-5. Стартовые структуры для этого шага являются выходами Шага 2.
4. Дальнейшее уточнение с использованием метода ЗМ с большим базовым набором
ПРИМЕЧАНИЕ: Цель здесь заключается в дальнейшем уточнении конфигурационной выборки кластеров Gly (H2O)n-1-5 с использованием лучшего квантово-механического описания. Стартовые структуры для этого шага являются выходами Шага 3.
5. Окончательные расчеты энергии и термодинамической коррекции
ПРИМЕЧАНИЕ: Цель здесь состоит в том, чтобы получить колебательное строение и энергии gly (H2O)n'1-5 кластеров с использованием большого базиса и ультратонкой интеграционной сетки для вычисления желаемых термохимических коррекций.
6. Вычисление атмосферных концентраций скоплений Gly (H2O)n'0-5 при комнатной температуре на уровне моря
ПРИМЕЧАНИЕ: Это достигается путем первого копирования термодинамических данных, полученных в предыдущем шаге в электронную таблицу и расчета Гиббс свободных энергий последовательной гидратации. Затем свободные энергии Гиббса используются для расчета равновесных констант для каждой последовательной гидратации. Наконец, набор линейных уравнений решается, чтобы получить равновесную концентрацию гидратов для данной концентрации мономеров, температуры и давления.
Первый набор результатов этого протокола должен быть набором низкоэнергетических структур Gly (H2O)n'1-5, найденных в рамках конфигурационной процедуры отбора проб. Эти структуры были оптимизированы на уровне теории PW91/6-311 и считаются точными для целей настоящего документа. Нет никаких доказательств того, что PW91/6-311»G» постоянно недооценивает или переоценивает связывающую энергию этих кластеров. Его способность предсказывать связывающие энергии относительно MP2/CBS32 и «DLPNO-CCSD»(T)/CBS60,61 оценки и эксперимент52 показывает много колебаний. То же самое относится и к большинству других функций плотности. Как правило, каждое значение n no 1 - 5 должно дать несколько низкоэнергетических структур в пределах около 5 ккал мол-1 из низкой энергетической структуры. Здесь мы сосредоточимся на первой структуре, произведенной скриптом run-thermo-pw91.csh для краткости. На рисунке 3 показаны самые низкие электронные энергетические изомеры кластеров Gly (H2O)n'0-5. Видно, что сеть водородных связей растет в сложности по мере увеличения числа молекул воды, и даже переходит от преимущественно планарной сети к трехмерной клетке- как структура на n No 5. В остальной части этого текста используются энергии и термодинамические количества, соответствующие этим пяти конкретным скоплениям.
Таблица 1 содержит термодинамические количества, необходимые для выполнения протокола. В таблице 2 показан пример выхода сценария run-thermo-pw91.csh, где печатаются электронные энергии, вибрационные коррекции нулевой точки и термодинамические коррекции при трех различных температурах. Для каждого кластера (строки) E'PW91/6-311-G, соответствует электронной энергии газовой фазы на уровне PW91/6-311'G, рассчитанной на сверхтонких интеграционных сетках в единицах Hartree, а также в нулевую точку вибрационной энергии(PVE) в единицах ккал mol-1. При каждой температуре, 216,65 K, 273.15 K, и 298,15 K, термодинамические коррекции перечислены, ЗХ enthalpy формирования в единицах ккал мол-1, S энтропии формирования в единицах cal mol-1, и Г Гиббс свободной энергии формирования в единицах кккал мол-1. В таблице 3 показанпримерный расчет общего количества свободного изменения энергии Гиббса гидратации, а также последовательного гидратации. Пример вычисления полного изменения энергии Gibbs свободной гидратации для реакции
начинается с расчета электронной энергии EPW91 как
где EPW91(H2O) берется из колонки Таблицы 2 C, а EPW91«Гли» и EPW91» H2O » взяты из столбца Таблицы 1 B. Далее мы вычисляем общее изменение энергетической точки газаE(0), включив изменение в нулевую точку колебательной энергии реакции как
получить колонку D. Здесь, «EPW91/6»311 »G» берется из столбца таблицы 3 C, E«PVE»Gly » (H2O)» из столбца D, и E«PVE»и E«PVE» H2O » из столбца таблицы 1 C. Ради краткости, мы перейдем к кластерам комнатной температуры, поэтому мы пропустим данные 216,65 K и 273.15 K. При комнатной температуре мы затем вычисляем enthalpy изменение реакцииH путем коррекции изменения энергии газовой фазы как
гдеиз столбца Таблицы 3 взята точка в год 3,h(H2O) берется из столбца таблицы 2 K, аH«Gly» и«H»H2O » взяты из столбца таблицы 1 J. Наконец, мы вычисляем Гиббса свободной изменения энергии реакции »G как
гдеиз столбца Таблицы 3 i, S«Gly» (H2O)) берется из столбца таблицы 2 L, а S«Gly» иSS2O » взяты из столбца таблицы 1 K. Обратите внимание, что значения энтропии должны быть преобразованы в единицы ккал-мола-1 K-1 во время этого шага. S
Теперь у нас есть необходимые количества для расчета атмосферных концентраций гидратированных глицинов, как показано в шаге 6. Результаты должны напоминать данные, приведенные в таблице 4,но следует ожидать небольших численных различий. В таблице 4 показаны равновесные концентрации гидратов, найденные в результате разработки системы шести уравнений в шаге 6.2 в одно матричного уравнения и его последующего решения. Мы начинаем с признания того факта, что система уравнений может быть написана как
где Kn является равновесной постоянной для nй последовательной гидратации глицина, w является концентрация воды в атмосфере, г является начальной концентрацией изолированных глицин в атмосфере, и gn является равновесной концентрацией Gly (H2O)n. Если мы переписываем вышеупомянутое уравнение как Ax q b,мы получим−1 х 1b, где ANo1 является обратной матрицы A. Это обратное можно легко вычислить с помощью встроенных функций электронной таблицы, как показано в таблице 4 для получения окончательных результатов.
На рисунке 4 показана равновесная концентрация гидратного глицина, рассчитанная в таблице 4 как функция температуры при 100% относительной влажности и 1 атмосферном давлении. Это показывает, что, как температура снижается с 298.15K до 216.65K, концентрация неивлажного глицина (n'0) уменьшается и у гидратированных глицина увеличивается. Глицин дигидрат (n'2), в частности, резко возрастает при снижении температуры, в то время как изменение концентрации других гидратов менее заметно. Эти обратные корреляции между температурой и концентрацией гидратов согласуется с ожиданием, что более низкие энергии Гиббса свободных гидратов при более низких температурах благоприятствуют образованию гидратов.
Рисунок 5 иллюстрирует относительную зависимость влажности равновесной концентрации гидратов глицина при 298,15К и 1 атмосферном давлении. Это ясно показывает, что по мере увеличения RH с 20% до 100%, концентрация гидратов (n'gt;0) увеличивается за счет негидратированных глицинов (n'0). Вновь прямая корреляция между относительной влажностью и концентрацией гидратов согласуется с мыслью о том, что наличие большего количества молекул воды при более высоком RH способствует образованию гидратов.
Как представлено, этот протокол дает качественное понимание гидратированных популяций глицина в атмосфере. Предполагая, что первоначальная концентрация изолированного глицина составляет 2,9 миллиона молекул на кубический сантиметр, мы видим, что негидратированный глицин (n'0) является наиболее распространенным видом в большинстве условий, за исключением T 216.65K и RH'100%. Дигидрат (n'2), который имеет самую низкую последовательную энергию гидратации Гиббса при всех трех температурах, является наиболее распространенным гидратом в условиях, рассматриваемых здесь. По прогнозам, моногидраты (n-1) и более крупные гидраты (n-3) будут найдены в незначительных количествах. При осмотре рисунка 3,обилие n и 1-4 кластеров может быть связано со стабильностью и напряжением в сети водородных связей кластеров. Эти скопления имеют молекулы воды водорода связаны с carboxylic кислоты moiety глицина в геометрии, тесно напоминающей те из различных водородных связанных кольцевых структур, что делает их особенно стабильными.
Рисунок 1: Схематическое описание текущей процедуры. Большой пул догадок, генерируемых генетическим алгоритмом (GA), уточняется с помощью серии оптимизаций геометрии PW91 до тех пор, пока не будет получен набор конвергентных структур. Вибрационные частоты этих структур вычисляются и используются для вычисления свободной энергии образования Гиббса, которая, в свою очередь, используется для расчета равновесных концентраций кластеров в условиях окружающей среды. Пожалуйста, нажмите здесь, чтобы просмотреть большую версию этой цифры.
Рисунок 2: Структура каталога представителей для каждого кластера. Внутренние скрипты, включенные в этот протокол, требуют приведенной выше структуры каталога, где n — это количество молекул воды. Для каждого n в gly-h2o-n, Есть следующие субдиректора: GA для генетического алгоритма с каталогом GA/pm7, «ММ для квантовой механики с «M/pw91-sb» для PW91/6-31-G, «M/pw91», «M/pw91»/ultrafine для оптимизации и окончательных вибрационных вычислений на сверхтонких сетках. Пожалуйста, нажмите здесь, чтобы просмотреть большую версию этой цифры.
Рисунок 3: Представитель низкоэнергетических структур Gly (H2O)n'0-5. Эти кластеры были электронной энергетической глобальной минимальной, оптимизированной на уровне теории PW91/6-311. Пожалуйста, нажмите здесь, чтобы просмотреть большую версию этой цифры.
Рисунок 4: Температурная зависимость Gly (H2O)n'0-5 как 100% относительная влажность и давление 1 атм. Концентрация гидратов дается в единицах молекул см-3. Пожалуйста, нажмите здесь, чтобы просмотреть большую версию этой цифры.
Рисунок 5: Относительная зависимость влажности Gly (H2O)n'0-5 как 298.15 K и давление 1 атм. Концентрация гидратов дается в единицах молекул см-3. Пожалуйста, нажмите здесь, чтобы просмотреть большую версию этой цифры.
Е-PW91/6-311 | 216.65 K | 273.15 K | 298.15 K | ||||||||
LB-UF | ЗПВЕ | ЗЧ | S | ЗГ | ЗЧ | S | ЗГ | ЗЧ | S | ЗГ | |
Воды | -76.430500 | 13.04 | 1.72 | 42.59 | 5.54 | 2.17 | 44.44 | 3.08 | 2.37 | 45.14 | 1.96 |
Глицин | -284.434838 | 48.55 | 2.65 | 69.53 | 36.14 | 3.70 | 73.81 | 32.09 | 4.22 | 75.61 | 30.22 |
Таблица 1: Мономерные энергии. Электронные энергии находятся в единицах Hartree в то время как все остальные количества находятся в единицах ккал мол-1. Вычислили воду и глицин на уровне теории и вибрационных частот PW91/6-311.G. С помощью thermo.pl скрипта были вычислены термодинамические коррекции давления 1 атм и температуры 298,15 К.
Е-PW91/6-311 | 0 K | 216.65 K | 273.15 K | 298.15 K | ||||||||
N | Имя | LB-UF | ЗПВЕ | ЗЧ | S | ЗГ | ЗЧ | S | ЗГ | ЗЧ | S | ЗГ |
1 | gly-h2o-1 | -360.88481 | 63.96 | 3.61 | 80.12 | 50.22 | 5.12 | 86.27 | 45.52 | 5.85 | 88.83 | 43.33 |
2 | gly-h2o-2 | -437.33763 | 79.33 | 4.53 | 90.86 | 64.17 | 6.46 | 98.78 | 58.81 | 7.40 | 102.06 | 56.30 |
3 | gly-h2o-3 | -513.78620 | 94.52 | 5.67 | 105.08 | 77.42 | 8.08 | 114.94 | 71.19 | 9.23 | 119.00 | 68.27 |
4 | gly-h2o-4 | -590.23667 | 109.80 | 6.03 | 104.98 | 91.30 | 8.78 | 116.21 | 84.40 | 10.11 | 120.87 | 81.14 |
5 | gly-h2o-5 | -666.68845 | 125.80 | 7.26 | 121.70 | 106.69 | 10.47 | 134.83 | 99.44 | 12.01 | 140.24 | 96.00 |
Таблица 2: Кластерные энергии. Энергии низшей энергии Gly (H2O)n'1-5 структур, найденных с помощью нашей процедуры, изложенной на рисунке 1. Электронные энергии находятся в единицах Hartree в то время как все остальные количества находятся в единицах ккал мол-1.
Полная гидратация: Гли nH2O йlt;--gt; Gly (H2O)n | Секвенциальная гидратация: Гли (H2O)n-1 - H2O lt;-gt; Gly (H2O)n | ||||||||||||||||
Е-PW91/6-311 | 216.65 | 273.15 | 298.15 | 216.65 | 273.15 | 298.15 | |||||||||||
N | системное имя | LB-UF | Зе (0) | ЗГ (T) | ЗГ (T) | ЗГ (T) | ЗГ (T) | ЗГ (T) | ЗГ (T) | LB-UF | Зе (0) | ЗГ (T) | ЗГ (T) | H (T) | ЗГ (T) | ЗГ (T) | ЗГ (T) |
1 | gly-h2o-1 | -12.22 | -9.85 | -10.61 | -3.68 | -10.61 | -1.87 | -10.59 | -1.07 | -12.22 | -9.85 | -10.61 | -3.68 | -10.61 | -1.87 | -10.59 | -1.07 |
2 | gly-h2o-2 | -26.22 | -21.53 | -23.10 | -9.27 | -23.11 | -5.66 | -23.09 | -4.06 | -14.00 | -11.68 | -12.49 | -5.59 | -12.50 | -3.79 | -12.50 | -2.99 |
3 | gly-h2o-3 | -37.56 | -30.72 | -32.88 | -12.90 | -32.87 | -7.69 | -32.82 | -5.38 | -11.34 | -9.19 | -9.78 | -3.63 | -9.76 | -2.03 | -9.73 | -1.32 |
4 | gly-h2o-4 | -50.10 | -40.34 | -43.48 | -15.87 | -43.54 | -8.71 | -43.51 | -5.55 | -12.54 | -9.62 | -10.60 | -2.97 | -10.67 | -1.02 | -10.69 | -0.17 |
5 | gly-h2o-5 | -63.45 | -51.41 | -55.42 | -20.58 | -55.51 | -11.48 | -55.48 | -7.45 | -13.35 | -11.07 | -11.94 | -4.71 | -11.97 | -2.77 | -11.97 | -1.90 |
Таблица 3: Энергии гидратации. Общая энергия гидратации и энергия последовательной гидратации для Gly (H2O)n'1-5 в единицах ккал мол-1. Здесь, E'PW91/6-311 "ГЗ" является изменение в электронной энергии, "E(0) является нулевой точки колебательной энергии (ЗПВЕ) исправленные изменения в энергии, H(T) является enthalpy изменения при температуре T, и "G"T) является Гиббс свободного изменения энергии гидратации каждого Gly (H2O)n-1.
Равномерное распределение гидратов в зависимости от температуры и относительной влажности | |||||||||
T 298.15K | T 273.15K | T 216.65K | |||||||
Гли (H2O)n | RH 100% | RH-50% | RH-20% | RH 100% | RH-50% | RH-20% | RH 100% | RH-50% | RH-20% |
0 | 1.3E-06 | 2.2E-06 | 2.7E-06 | 1.1E-06 | 2.0E-06 | 2.7E-06 | 6.1E-05 | 1.5E-06 | 2.5E-06 |
1 | 2.3E-05 | 1.9E-05 | 9.5E-04 | 2.0E-05 | 1.9E-05 | 9.9E-04 | 1.2E-05 | 1.5E-05 | 9.5E-04 |
2 | 1.0E-06 | 4.3E-05 | 8.4E-04 | 1.3E-06 | 6.1E-05 | 1.3E-05 | 1.8E-06 | 1.1E-06 | 3.0E-05 |
3 | 2.8E-05 | 5.8E-04 | 4.5E-03 | 3.2E-05 | 7.4E-04 | 6.3E-03 | 3.1E-05 | 9.6E-04 | 1.0E-04 |
4 | 1.1E-04 | 1.1E-03 | 3.4E-01 | 1.3E-04 | 1.5E-03 | 5.0E-01 | 1.1E-04 | 1.8E-03 | 7.5E-01 |
5 | 7.5E-03 | 3.9E-02 | 4.9E-00 | 1.2E-04 | 7.2E-02 | 9.7E-00 | 2.4E-04 | 1.9E-03 | 3.1E-01 |
Таблица 4: Равновесные увлажняемые концентрации Gly (H2O)n'0-5 в качестве функциональной температуры (T'298.15K, 273.15K, 216.65K) и относительная влажность (RH-100%, 50%, 20%). Концентрация гидратов дается в единицах молекул см-3 при условии экспериментальных значений56,,57,,58, из «Gly»0 , 2,9 х 106 см-3 и Х2O » 7,7 х 1017 см-3, 1,6 х 1017 см-3 и 9,9 х 1010 см-3 при 100% относительной влажности, 273.15 K и 216.65 K, соответственно59.
Дополнительные файлы. Пожалуйста, нажмите здесь, чтобы загрузить эти файлы.
Точность данных, генерируемых этим протоколом, зависит главным образом от трех моментов: i) разнообразие конфигураций, отобранных шагом 2, (ii) точность юндаром структуры системы (iii) и точность термодинамических коррекций. Каждый из этих факторов может быть решен путем изменения метода путем редактирования включенных скриптов. Первый фактор легко преодолеть с помощью большего первоначального пула случайно генерируемых структур, более многочисленных итераций ГА и более слабого определения критериев, участвующих в ГА. Кроме того, можно использовать другой полуэмпирический метод, такой как самосогласованная плотность заряда-функциональная герметичная герметичная (SCC-DFTB)62 модель и эффективный фрагмент потенциал (EFP)63 модель для того, чтобы исследовать эффекты различных физических описаний. Основным ограничением здесь является неспособность метода формировать или ломать ковалентные связи, а это означает, что мономера заморожены. Процедура ГА находит только наиболее стабильные относительные позиции этих замороженных мономеров в соответствии с полуэмпирическим описанием.
Точность электронной структуры системы может быть улучшена различными способами, каждый из которых имеет свои вычислительные затраты. Можно выбрать более высокую плотность функциональных, таких как M06-2X64 и wB97X-V65, или квантово-механических (ЗМ) метод, таких как Мюллер-Плессет66,67,68 (MPn) возмущения теорий и совместно-кластер69 (CC) методы для улучшения физического описания системы. В иерархии функциональных функций, производительность, как правило, улучшается при переходе от обобщенного градиента приближения (GGA) функций, как PW91 в диапазоне разделенных гибридных функций, как wB97X-D и мета-GGA гибридных функций, как M06-2X.
Недостатком методов DFT является то, что систематическое сближение с точным значением невозможно; однако, методы DFT являются вычислительно недорогими и существует широкий спектр функций для широкого спектра приложений.
Энергия, рассчитанная с использованием методов волнообразов, таких как MP2 и CCSD(T) в сочетании с корреляционными последовательными базовыми наборами увеличения кардинального числа («Ауг-кк-пВ»,Т,...» В) сходятся к их полной основе установленного предела систематически, но вычислительная стоимость каждого расчета становится непомерно высокой по мере роста размера системы. Дальнейшее уточнение электронной структуры может быть достигнуто с помощью явно коррелированных базовых наборов70 и экстраполации на полный базовый набор (CBS)71 предел. Наша недавняя работа показывает, что плотность установлены явно коррелирует второго порядка Мюллер-Плессет (DF-MP2-F12) возмущенный подход дает энергии приближается к MP2 / CBS вычислений32. Изменение текущего протокола для использования различных методов электронной структуры включает в себя два шага: (i) подготовить шаблон входиного файла после синтаксиса, приведенного программным обеспечением, (ii) и откорректировать run-pw91-sb.csh, run-pw91-lb.csh, и run-pw91-lb-ultrafine.csh скрипты для создания правильного ввода.
Наконец, точность термодинамических коррекций зависит от метода электронной структуры, а также от описания PES вокруг глобального минимума. Точное описание PES требует вычисления производных ПЭУ третьего и более высокого порядка в отношении перемещений в ядерных степенях свободы, таких как поле квартической силы72,,73 (КФФ), что является исключительно дорогостоящей задачей. Текущий протокол использует гармоничное приближение осциллятора к вибрационным частотам, в результате чего необходимо вычислить только до вторых производных PES. Этот подход становится проблематичным в системах с высокой агармоничностью, таких как очень дискеты молекул и симметричных двойного потенциала из-за большой разницы в истинном PES и гармонических PES. Кроме того, стоимость высококачественного ПЭУ от вычислительно-требовательного метода электронной структуры только усугубляет проблему затрат на вибрационные частотные вычисления. Один из подходов к преодолению этого заключается в использовании электронных энергий от высококачественного электронного расчета структуры наряду с вибрационными частотами, вычисляемыми на более низком качестве PES, что приводит к балансу между стоимостью и точностью. Текущий протокол может быть изменен с целью использования различных описаний ПЭУ, описанных в предыдущем пункте; однако, можно также отсеивать биевания частотных ключевых слов в сценариях и шаблонах для вычисления агармонических вибрационных частот.
Двумя важнейшими вопросами для любого протокола выборки конфигурации являются первоначальный метод отбора проб потенциальной энергетической поверхности и критерии, используемые для определения каждого кластера. Мы широко использовали различные методы в нашей предыдущей работе. Для первого выпуска, первоначального метода отбора проб потенциальной энергетической поверхности, мы сделали выбор использования ГА с полуэмпирическими методами на основе этих факторов. Конфигурационная выборка с использованием химической интуиции26, случайная выборка, и молекулярной динамики (MD)29,30, не в состоянии найти ориентировочный глобальной минимы регулярно для кластеров больше, чем 10 мономеров, как мы наблюдали в наших исследованиях водных кластеров18. Мы успешно использовали бассейн прыжков (BH) для изучения сложных PES (H2O)1174, но это требует ручного включения некоторых потенциальных изомерителей низкой энергии алгоритм BH не нашли. Сравнение производительности BH и GA в поиске глобального минимума водных кластеров, (H2O)n'10-20 показали, что GA последовательно нашел глобальный минимум быстрее, чем BH75. GA, как реализовано в OGOLEM и CLUSTER является очень универсальным, потому что он может быть применен к любому молекулярному кластеру, и он может взаимодействовать с огромным количеством пакетов с классическим силовым полем, полуэмпирическими, функциональными плотностями и возможностями ab initio. Выбор PM7 определяется его скоростью и разумной точностью. Практически любой другой полуэмпирический метод будет иметь значительно более высокую вычислительную стоимость.
Что касается второго выпуска, мы изучили с использованием различных критериев для выявления уникальных структур, начиная от электронных энергий, дипольные моменты, перекрытия RMSDs и вращательных констант. Использование дипольные моменты оказалось трудным, потому что оба компонента диполи-момента зависели от ориентации молекулы, а общий дипольные моменты были очень чувствительны к геометрии различий таким образом, что было трудно установить пороги, определяющие структуры, одинаковые или уникальные. Сочетание электронных энергий и вращательных констант оказалось наиболее полезным.
Текущие критерии для того, чтобы считать две структуры уникальными, основаны на пороге разницы в энергии 0,10 ккалмол -1 и постоянной разнице в повороте в 1%. Таким образом, две структуры считаются различными, если их энергия отличается более чем на 0,10 ккал мол-1 (0,00015 a.u.) И любая из трех их вращательных констант (A, B, C) отличается более чем на 1%. Существенные внутренние ориентиры на протяжении многих лет сочли эти пороговые значения разумным выбором. Наш конфигурационный подход к выборке и методология скрининга были применены к очень слабо связанным кластерам, таким как полиароматические углеводороды, комплексированные с водой76,,77, а также сильно связанные тернарисные гидраты сульфата, содержащие аммиак и амины32. Для кластеров, где существуют различные состояния протонации, который следует учитывать, лучшим подходом является запуск различных расчетов ГА, каждый из которых начинается с мономеров в разных состояниях протонации. Это гарантирует тщательное обдумание структур с различными состояниями протонации. Однако низкоуровневые расчеты DFT часто позволяют состояниям протонации изменяться в ходе оптимизации геометрии, тем самым уступая наиболее стабильное состояние протонации независимо от стартовой геометрии.
Наши методы гамма-конфигурации выборки должны хорошо работать даже для дискетных молекул, пока коды GA взаимосвязаны с общими, непараметризированными методами, которые позволяют мономерам принимать различные конфигурации в ходе бега ГА. Например, взаимодействие ГА с PM7 позволит изменить структуры мономеров, но если их облигации прервутся, как это произойдет при смене состояний протонации, структуры могут быть отвергнуты как неприемлемые кандидаты.
Мы рассмотрели различные способы исправления недостатков гармонического приближения, особенно те, которые возникают из-за низких колебатевых частот. Включение квазигармонического приближения в нынешнюю методологию не сложно. Тем не менее, есть еще вопросы о квази-гармонический метод, особенно когда дело доходит до частоты отсечения, ниже которого он будет применяться. Кроме того, нет строгих контрольных работ, изучающих надежность приближения квази-RRHO, хотя общепринятая точка зрения предполагает, что это должно быть улучшением по сравнению с приближением RRHO.
Представленный таким образом протокол может быть обобщен любой системе молекулярных кластеров нековалентно связанных газовых фаз. Он также может быть обобщен использовать любой полуэмпирический метод, электронный метод структуры и программного обеспечения, а также метод вибрационного анализа и программного обеспечения путем редактирования скриптов и шаблонов. Это предполагает, что пользователь доволен интерфейсом командной строки Linux, скриптом Python и высокопроизводительными вычислениями. Незнакомый синтаксис и внешний вид операционной системы Linux и отсутствие опыта работы с криптами является самой большой ловушкой в этом протоколе и где новые студенты борются больше всего. Этот протокол успешно используется в различных реализациях в течение многих лет в нашей группе, в основном упором на влияние серной кислоты и аммиака на образование аэрозолей. Дальнейшие усовершенствования этого протокола будут включать в себя более надежный интерфейс для более электронной структуры программного обеспечения, альтернативные реализации генетического алгоритма, и, возможно, использование новых методов для более быстрых вычислений электронных и вибрационных энергий. Наши текущие применения этого протокола изучают важность аминокислот на ранних стадиях формирования аэрозолей в текущей атмосфере и в формировании более крупных биологических молекул в пребиотической среде.
Ни один.
Этот проект был поддержан грантами CHE-1229354, CHE-1662030, CHE-1721511 и CHE-1903871 от Национального научного фонда (GCS), Арнольда и Мэйбл Бекман Фонд Бекман Стипендиат премии (AGG), и Барри М. Голдуотер стипендии (AGG). Использовались высокопроизводительные вычислительные ресурсы консорциума MERCURY (http://www.mercuryconsortium.org).
Name | Company | Catalog Number | Comments |
Avogadro | https://avogadro.cc | Open-source molecular visualization program | |
Gaussian [09/16] Software | http://www.gaussian.com/ | Commercial ab initio electronic structure program | |
MOPAC 2016 | http://openmopac.net/MOPAC2016.html | Open-source semi-empirical program | |
OGOLEM Software | https://www.ogolem.org | Genetic algorithm-based global optimization program | |
OpenBabel | http://openbabel.org/wiki/Main_Page | Open-source cheminformatics library | |
calcRotConsts.py | Shields Group, Department of Chemistry, Furman University | Python script to compute rotational constants | |
calcSymmetry.csh | Shields Group, Department of Chemistry, Furman University | Shell script to calculate symmetry number of a molecule given Cartesian coordinates | |
combine-GA.csh | Shields Group, Department of Chemistry, Furman University | Shell script to combine energy and rotational constants from different GA directories | |
combine-QM.csh | Shields Group, Department of Chemistry, Furman University | Shell script to combine energy and rotational constants from different QM directories | |
gaussianE.csh | Shields Group, Department of Chemistry, Furman University | Shell script to extract Gaussian 09 energies | |
gaussianFreqs.csh | Shields Group, Department of Chemistry, Furman University | Shell script to extract Gaussian 09 vibrational frequencies | |
getrotconsts | Shields Group, Department of Chemistry, Furman University | Executable to calculate rotational constants given a molecule's Cartesian coordinates | |
getRotConsts-dft-lb.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute rotational constants for a batch of large basis DFT optimized structures | |
getRotConsts-dft-lb-ultrafine.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute rotational constants for a batch of ultrafine DFT optimized structures | |
getRotConsts-dft-sb.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute rotational constants for a batch of small basis DFT optimized structures | |
getRotConsts-GA.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute rotational constants for a batch of genetic algorithm optimized structures | |
global-minimum-coords.xyz | Shields Group, Department of Chemistry, Furman University | Cartesian coordinates of global minimum structures of gly-(h2o)n, where n=0-5 | |
make-thermo-gaussian.csh | Shields Group, Department of Chemistry, Furman University | Shell script to extract data from Gaussian output files and make input files for the thermo.pl script | |
ogolem-input-file.ogo | Shields Group, Department of Chemistry, Furman University | Ogolem sample input file | |
ogolem-submit-script.pbs | Shields Group, Department of Chemistry, Furman University | PBS batch submission file for Ogolem calculations | |
README.docx | Shields Group, Department of Chemistry, Furman University | Clarifications to help readers use the scripts effectively | |
runogolem.csh | Shields Group, Department of Chemistry, Furman University | Shell script to run OGOLEM | |
run-pw91-lb.csh | Shields Group, Department of Chemistry, Furman University | Shell script to run a batch of large basis DFT optimization calculations | |
run-pw91-lb-ultrafine.csh | Shields Group, Department of Chemistry, Furman University | Shell script to run a batch of ultrafine DFT optimization calculations | |
run-pw91-sb.csh | Shields Group, Department of Chemistry, Furman University | Shell script to run a batch of small basis DFT optimization calculations | |
run-thermo-pw91.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute the thermodynamic corrections for a batch of DFT optimized structures | |
similarityAnalysis.py | Shields Group, Department of Chemistry, Furman University | Python script to determine unique structures based on rotational constants and energies | |
symmetry | Shields Group, Department of Chemistry, Furman University | Executable to calculate molecular symmetry given Cartesian coordinates | |
symmetry.c | (C) 1996, 2003 S. Patchkovskii, Serguei.Patchkovskii@sympatico.ca | C code to determine the molecular symmstry of a molecule given Cartesian coordinates | |
template-marcy.pbs | Shields Group, Department of Chemistry, Furman University | Template for a PBS submit script which uses OGOLEM | |
template-pw91.com | Shields Group, Department of Chemistry, Furman University | Template Gaussian 09 input | |
template-pw91-HL.com | Shields Group, Department of Chemistry, Furman University | Template Gaussian 09 input for ultrafine DFT optimization | |
thermo.pl | https://www.nist.gov/mml/csd/chemical-informatics-research-group/products-and-services/program-computing-ideal-gas | Perl open-source script to compute ideal gas thermodynamic corrections | |
gly-h2o-n.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet for the complete protocol | |
table-1.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet | |
table-2.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet | |
table-3.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet | |
table-4.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet | |
water.xyz | Shields Group, Department of Chemistry, Furman University | Cartesian coordinates of water | |
glycine.xyz | Shields Group, Department of Chemistry, Furman University | Cartesian coordinates of glycine |
Запросить разрешение на использование текста или рисунков этого JoVE статьи
Запросить разрешениеThis article has been published
Video Coming Soon
Авторские права © 2025 MyJoVE Corporation. Все права защищены