2023
2022
2021
2020
2019
2018
2017
2016
Квантово-химические рассчеты для построения топологии: гессиан во внутренних координатах
В этой задаче мы должны будем получить гессиан во внутренних координатах. Внутренние координаты молекулы - некоторый набор невырожденных координат, который кажется исследователю более естественным, чем декартовы координаты. Например, всегда можно построить полный координатный базис только из длин связей, величин валентных и торсионных углов. Диагональные элементы матрицы Гессе, записанной в такой системе координат имеют смысл жесткости соответствующей связи или угла (с некоторыми оговорками, происходящими из-за неоднозначного набора углов).
Итак, на сервере kodomo установлена программа FireFly Александра Грановского, которая является ответлением пакета US-GAMESS. От него она наследует весь синтаксис. В тех местах, где Firefly сильно отличается от US-GAMESS изменения в синтаксисе значительны. Сейчас программа установлена в /Z/local/ff81, а ее запуск производится запускающим скриптом:
1 /Z/local/bin/pcgamess2.py -n 2 -i input.inp -o OUT.log &
Где n - количество ядер, на которых FireFly следует распараллелить.
I. На первом этапе вы должны построить молекулу и оптимизировать ее геометрию.
Молекулу вы выбираете и строите самостоятельно теми средствами, которыми хотите. Молекула не должна быть очень большой. Следует ограничиться парой десятков атомов. Координаты для FireFly вы можете получить babel-ем, например:
1 babel -ipdb file.pdb -ogamin file.inpСконструируйте свой INPUT на основе следующего шаблона. Обратите внимание на заряд, мультиплетность, тип волновой функции (SCFTYP) и наличие диффузных функций в базисе.
1 $contrl 2 runtyp=optimize 3 ! RHF if mult=1, UHF or ROHF if mult>1 4 scftyp=RHF 5 DFTTYP=B3LYP5 6 ! your charge and multiplicity 7 icharg=0 mult=2 8 $end 9 $basis gbasis=n31 ngauss=6 npfunc=1 ndfunc=1 10 ! use diffuse functions (.t.) for anions! false (.f.) in other case. 11 diffsp=.f. 12 $end 13 $SCF FDIFF=.f. DIIS=.t. DIRSCF=.t. $END 14 ! always use following groups for smp machine 15 $p2p p2p=.t. dlb=.t. $end 16 $smp csmtx=.t. call64=.t. $end
Проведите расчет. Проверьте, что он завершился корректно. Посмотрите на то, как происходила минимизация энергии, используя программу GabEdit. Проверьте, что расчет сошелся - тогда в output-файле присутствует строка "EQUILIBRIUM GEOMETRY LOCATED". Если такой строки нет - переходите к №4.
Для продолжения расчета вы можете забрать последние координаты через babel -igamout, а можете прямо из output-файла, в месте последнего NSERCH-а:
1 ... 2 NSERCH= 5 3 4 5 COORDINATES OF ALL ATOMS ARE (ANGS) 6 ATOM CHARGE X Y Z 7 ------------------------------------------------------------ 8 C 6.0 0.7749819844 -0.0036804832 -0.0000311613 9 C 6.0 -0.7159679493 0.0003550144 0.0000301921 10 H 1.0 1.3296953993 0.9293635358 -0.0002162206 11 H 1.0 1.3428886995 -0.9273597055 -0.0001174462 12 H 1.0 -1.1283971469 -1.0132774776 -0.0021729884 13 H 1.0 -1.1234133599 0.5234315680 0.8794327568 14 H 1.0 -1.1234436209 0.5274378673 -0.8769145615 15 ...
Эти координаты могут быть напрямую вставлены в блок $DATA.
II. На втором этапе нам нужно провести более глубокую минимизацию геометрии для того, чтобы в гессиане не было мнимых частот. Глубокую оптимизацию лучше провести во внутренних координатах.
Постройте Z-матрицу в формате gaussian при помощи babel -ogzmat, используя в качестве координат финальные координаты минимизации геометрии на предыдущем этапе. Их можно забрать babel-ем из output-файла FireFly с флагом -igamout. Полученную матрицу вставьте в блок $DATA и составьте check-расчет ("EXETYP=CHECK" в группе $CONTRL). Полученный output-файл сконвертируйте во что-нибудь и посмотрите глазами. Блок $data должен получиться похожим на это:
1 $data 2 Ethyl radical 3 C1 4 C 5 C 1 r2 6 H 1 r3 2 a3 7 H 1 r4 2 a4 3 d4 8 H 2 r5 1 a5 3 d5 9 H 2 r6 1 a6 3 d6 10 H 2 r7 1 a7 3 d7 11 12 r2 1.4587 13 r3 1.0814 14 a3 120.92 15 r4 1.0807 16 a4 121.99 17 d4 180.03 18 r5 1.0971 19 a5 112.83 20 d5 179.87 21 r6 1.0976 22 a6 110.22 23 d6 300.35 24 r7 1.0976 25 a7 110.26 26 d7 59.36 27 $end
Для того, чтобы сам расчет происходил во внутренних координатах, необходимо составить внутренние координаты на базе Z-матрицы. Можно написать для этого пару строк на питоне, можно сделать вручную, можно в экселе. В итоге должен получиться код такого вида:
Первая цифра означает тип внутренней координаты (1 -связь, 2 - угол, 3 - двугранный угол), следующие цифры - номера атомов. Если группа $zmat задана, координаты могут быть указаны как угодно, в том числе и в декартовой системе. После оформления $zmat не забудьте установить "nzvar=<3N-6>" в группе $control.
Теперь ваш input-файл готов для проведения оптимизации во внутренних координатах. Запустите ее. Проверьте, что оптимизация действительно идет во внутренних координатах. В этом случае после каждого шага NSERCH в выходном файле должна печататься таблица внутренних координат:
1 -------------------- 2 INTERNAL COORDINATES 3 -------------------- 4 5 - - ATOMS - - COORDINATE COORDINATE 6 NO. TYPE I J K L M N (BOHR,RAD) (ANG,DEG) 7 ------------------------------------------------------------------------------ 8 1 STRETCH 2 1 2.7565433 1.4587000 9 2 STRETCH 3 1 2.0435497 1.0814000 10 3 STRETCH 4 1 2.0422269 1.0807000 11 4 STRETCH 5 2 2.0732184 1.0971000 12 5 STRETCH 6 2 2.0741632 1.0976000 13 6 STRETCH 7 2 2.0741632 1.0976000 14 7 BEND 3 1 2 2.1104521 120.9200000 15 8 BEND 4 1 2 2.1291272 121.9900000 16 9 BEND 5 2 1 1.9692550 112.8300000 17 10 BEND 6 2 1 1.9237019 110.2200000 18 11 BEND 7 2 1 1.9244000 110.2600000 19 12 TORSION 4 1 2 3 -3.1410691 -179.9700000 20 13 TORSION 5 2 1 3 3.1393237 179.8700000 21 14 TORSION 6 2 1 3 -1.0410889 -59.6500000 22 15 TORSION 7 2 1 3 1.0360274 59.3600000
- Проверьте, что ваша оптимизация прошла успешно. Изучите выходной файл, посмотрите на молекулу глазами.
III. Заключительный этап - нахождение гессиана во внутренних координатах.
Возьмите последние координаты последней оптимизации. Достаточно декартовых. Весь остальной инпут сделайте таким же, поставьте "runtype=hessian". Для перевода гессиана во внутренние координаты добавьте группу "$force prtifc=.t. $end". Запустите расчет.
Внимательно изучите выходной файл. Вы должны увидеть строки:
.. за которыми последует матрица. Кроме того, вы должны увидеть частоты, их приведенные массы и интенсивности - характеристики IR спектра.
Буква "I" после частоты 104.70 означает мнимую единицу. Эта колебательная мода - мнимая. В таком случае говорят о неистинно стационарном состоянии. Требуется более глубокая минимизация. В этом случае переходите к шагу 3, если все в порядке - к шагу 4.
Для более глубокой минимизации установите в групе $statpt параметры "opttol=1e-6" порог, до которого спускаться, и "nstep=200" - максимальное количество шагов спуска.
- Откройте програмой GABEDIT ваш гессиановый OUTPUT-файл, визуализируйте самую интенсивную нормальную моду, постройте ИК-спектр.