Entry tags:

cp2k, it, linux, лабораторія, наука

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

Все файлы я выкладываю на облачное хранилище MEGA, они будут доступны по ссылкам в постах. Если удобнее куда-то в другое место --- напишите.

В частности, файлы для данного расчета лежат здесь: https://mega.nz/#F!oQJ00bjY!wAFcNt6f5eCJIKFJczFFDg

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

Конфигурационный файл CP2K состоит из секций (вкладываемых друг в друга) и параметров. Начало секции объявляется амперсандром & с названием секции, завершение -- амперсандом со словом end и необязательным названием секции, например:
&global

&end global
Или:
&global

&end

Несмотря на то, что название секции после &end необязательно, его лучше указывать, чтобы не запутаться.

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

&global
print_level medium
&end global

Восклицательный знак используется для объявления комментариев (все содержимое от восклицательного знака до конца строки игнорируется). Регистр букв не имеет значения --- по традиции почему-то обычно пишут прописными, но я лично предпочитаю строчные.

Набор возможных настроек обширный, но настраивать, разумеется, нужно не все. Большинство параметров имеют относительно разумные значения по умолчанию. Таким образом, мы изменяем параметры там, где:
--- Мы непосредственно хотим поменять поведение программы;
--- Где значения по умолчанию нас не устраивают;
--- Где значения по умолчанию не совместимы с другими настройками расчета, которые нам нужны (к сожалению, CP2K не умеет менять одни параметры в зависимости от значений других).

Обычно для учебных расчетов пользователи CP2K почему-то предпочитают кластеры воды, но оставим это уныние и попробуем оптимизировать что-то повеселее --- например, молекулу LSD-25.

Для получения начальных координат молекулы воспользуемся свободной программой-редактором Avogadro, которую можно скачать с официального сайта или установить из репозитория. В Avogadro нужно нарисовать молекулу по структурной формуле, после чего оптимизировать ее геометрию молекулярной механикой (Меню -- Расширения -- Оптимизировать геометрические параметры, повторять пока молекула не перестанет изменяться), результат сохранить в файл в формате XYZ (*.xyz) (в каталоге с файлами расчета результат сохранен как LSD_mm.xyz).

Ниже привожу комментированный файл конфигурации (input-файл). Звездочки означают пояснения после текста файла.

&global                         ! Секция, отвечающая за глобальные настройки расчета
    project opt_LSD                     ! Основа для имени служебных файлов расчета
    run_type geo_opt                    ! Тип расчета -- оптимизация геометрии
&end global
&force_eval
    method quickstep                    ! Квантовохимический расчет
    &dft                                ! Секция с параметрами КМ-расчета, называется dft независимо от того, используется ли именно метод DFT
        &qs                                     ! Выбор квантовохимического метода
            method PM6                                  ! Используем полуэмпирический метод PM6
        &end qs
        &scf                                    ! Настройки процедуры итеративного решения уравнений Хартри-Фока или Кона-Шэма
            max_scf 50                                  ! Не более 50 шагов*
            &ot                                         ! Использовать алгоритм "orbital transformation"*
                preconditioner full_single_inverse              ! Предобуславливатель для алгоритма сопряженных градиентов*
            &end ot
            &outer_scf                                  ! Настройки "внешнего цикла" процедуры итеративного решения*
                max_scf 50                                      ! Не более 50 шагов*
            &end outer_scf
        &end scf
    &end dft
    &subsys                             ! Описание исследуемой системы
        &cell                                   ! Настройки ячейки
            ABC 20 20 20                                ! Размер ячейки (a,b,c)*
            periodic none                               ! Непериодическая система
        &end cell
        &topology
            &center_coordinates .TRUE.                  ! Центрировать молекулу в ячейке*
            &end center_coordinates
        &end topology
        &coord                                  ! Координаты атомов*
C         -4.44076        3.38549       -0.07518
C         -4.20245        2.00598       -0.18155
C         -3.40353        4.29868        0.16462
C         -2.08153        3.85820        0.30373
C         -1.85203        2.48061        0.19671
C         -2.87974        1.58242       -0.03758
H         -5.45878        3.75576       -0.18081
H         -5.01183        1.30889       -0.36837
H         -3.64539        5.35408        0.24965
N         -2.29621        0.33449       -0.08450
C         -0.93528        0.43399        0.11257
C         -0.63940        1.76582        0.28801
H         -0.30972       -0.44833        0.09907
C         -0.86508        4.65877        0.57755
C          0.36039        3.90614        1.19518
C          0.63147        2.50878        0.51398
H          1.29060        1.89916        1.14186
H          1.12929        2.65644       -0.45309
C         -0.82945        5.97524        0.29695
C          0.36951        6.81146        0.65051
C          1.34889        6.11639        1.60487
N          1.61372        4.72337        1.19341
H          0.99927        6.13518        2.64110
H          2.28636        6.68867        1.59278
H          0.07569        3.71300        2.24199
H         -1.65610        6.48525       -0.18794
C          2.65008        4.14824        2.05441
H          2.30988        4.03449        3.08985
H          3.54377        4.78317        2.05762
H          2.98291        3.17369        1.68576
C         -0.00730        8.22627        1.12929
H          0.90263        6.96361       -0.30101
O          0.34565        9.18339        0.42712
N         -0.72138        8.44549        2.30864
C         -1.08197        9.83731        2.60724
C         -2.35726       10.26645        1.89496
H         -0.26177       10.49877        2.30508
H         -1.19818        9.96276        3.68753
H         -2.25710       10.17414        0.80878
H         -2.58814       11.31039        2.12863
H         -3.20913        9.65264        2.20403
C         -1.33483        7.41700        3.15831
C         -0.76228        7.39834        4.57002
H         -1.21067        6.42494        2.73341
H         -2.41362        7.60479        3.19401
H          0.32770        7.30253        4.55621
H         -1.17476        6.55406        5.13152
H         -1.00874        8.31092        5.12024
H         -2.79741       -0.52733       -0.24867
        &end coord
    &end subsys
&end force_eval

Пояснения:
1. Для поиска решения систем уравнений Хартри-Фока или Кона-Шэма (решения самосогласованного поля) CP2K, в отличие от многих других программ, в большинстве случаев рекомендует отдельный алгоритм, называемый orbital transformation, хотя умеет пользоваться и "стандартными" алгоритмами, связанными с диагонализацией фокиана, в т.ч. DIIS. Алгоритм основан на прямом поиске локального минимума электронной энергии как функции орбитальных коэффициентов, рекомендуемым методом минимизации при этом является метод сопряженных градиентов с предобуславливателем. В принципе, тщательная настройка секции &ot нужна только в случае проблем со сходимостью решения самосогласованного поля, однако в данном случае мы вынуждены переопределить предобуславливатель, потому что предобуславливатель по умолчанию, full_kinetic, не реализован для полуэмпирических методов.
2. CP2K организует решение ССП (самосогласованного поля, или self-consistent field, SCF) с помощью т.н. внешних и внутренних циклов. Внутренний цикл -- непосредственно итерация решателя ССП. Внешний цикл служит для перезапуска внутренних (возможна более тонкая настройка, но в простейшем случае просто происходит перезапуск). Смысл этой акции в том, что при перезапуске внутренних циклов пересчитывается предобуславливатель, что существенно улучшает сходимость. Таким образом, в данном расчете может быть выполнено до 50 внешних циклов, каждый из которых может быть до 50 итераций. Реально ССП сойдется гораздо раньше, разумеется.

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

3. Ячейка и геометрия. Архитектура CP2K такова, что даже для непериодических расчетов необходимо задавать ячейку. Причем это не формальное требование --- исследуемая молекула должна помещаться в ячейке (иначе программа вылетит с ошибкой). Кроме того, подразумевается, что на границе ячейки электронная плотность будет нулевой (т.е. ниже порога отсекания) --- ошибки в противном случае не будет, но будет много злобных предупреждений UPD: верно только для вэйвлетного решателя уравнения Пуассона, в общем случае необязательно. В то же время, чрезмерное увеличение ячейки нежелательно, т.к. с ростом ячейки растут требования к вычислительным ресурсам. Поэтому задаем достаточно большую ячейку 20x20x20 ангстрем и (UPD: опционально) просим программу центрировать молекулу в этой ячейке. Секция &coords, задающая геометрию --- это список атомов с указанием элемента и декартовых координат в ангстремах (фактически, содержимое файла LSD_mm.xyz без первых двух строк --- числа атомов и комментария --- именно таким образом и следует заполнить эту секцию).

Итак, сохраняем input как opt_LSD.inp и запускаем:
cp2k.ssmp opt_LSD.inp | tee opt_LSD.out

Разбор полетов --- в следующих сериях.