Kodomo

Пользователь

Учебная страница курса биоинформатики,
год поступления 2010

Квантово-химические рассчеты для построения топологии: гессиан во внутренних координатах

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

Итак, на сервере 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. На первом этапе вы должны построить молекулу и оптимизировать ее геометрию.

  1. Молекулу вы выбираете и строите самостоятельно теми средствами, которыми хотите. Молекула не должна быть очень большой. Следует ограничиться парой десятков атомов. Координаты для FireFly вы можете получить babel-ем, например:

       1 babel -ipdb file.pdb -ogamin file.inp
    
  2. Сконструируйте свой 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
    
  3. Проведите расчет. Проверьте, что он завершился корректно. Посмотрите на то, как происходила минимизация энергии, используя программу GabEdit. Проверьте, что расчет сошелся - тогда в output-файле присутствует строка "EQUILIBRIUM GEOMETRY LOCATED". Если такой строки нет - переходите к №4.

  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. На втором этапе нам нужно провести более глубокую минимизацию геометрии для того, чтобы в гессиане не было мнимых частот. Глубокую оптимизацию лучше провести во внутренних координатах.

  1. Постройте 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
    
  2. Для того, чтобы сам расчет происходил во внутренних координатах, необходимо составить внутренние координаты на базе Z-матрицы. Можно написать для этого пару строк на питоне, можно сделать вручную, можно в экселе. В итоге должен получиться код такого вида:

       1  $zmat izmat(1)=
       2     1, 2, 1
       3     1, 3, 1,
       4     1, 4, 1,
       5     1, 5, 2,
       6     1, 6, 2,
       7     1, 7, 2,
       8     2, 3, 1, 2,
       9     2, 4, 1, 2,
      10     2, 5, 2, 1,
      11     2, 6, 2, 1,
      12     2, 7, 2, 1,
      13     3, 4, 1, 2, 3,
      14     3, 5, 2, 1, 3,
      15     3, 6, 2, 1, 3,
      16     3, 7, 2, 1, 3
      17  $end
    

    Первая цифра означает тип внутренней координаты (1 -связь, 2 - угол, 3 - двугранный угол), следующие цифры - номера атомов. Если группа $zmat задана, координаты могут быть указаны как угодно, в том числе и в декартовой системе. После оформления $zmat не забудьте установить "nzvar=<3N-6>" в группе $control.

  3. Теперь ваш 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
    
  4. Проверьте, что ваша оптимизация прошла успешно. Изучите выходной файл, посмотрите на молекулу глазами.


III. Заключительный этап - нахождение гессиана во внутренних координатах.

  1. Возьмите последние координаты последней оптимизации. Достаточно декартовых. Весь остальной инпут сделайте таким же, поставьте "runtype=hessian". Для перевода гессиана во внутренние координаты добавьте группу "$force prtifc=.t. $end". Запустите расчет.

  2. Внимательно изучите выходной файл. Вы должны увидеть строки:

       1  HESSIAN MATRIX IN INTERNAL COORDINATES
       2  UNITS ARE HARTREE/BOHR**2, HARTREE/RADIAN**2, OR HARTREE/BOHR*RADIAN
    
  3. .. за которыми последует матрица. Кроме того, вы должны увидеть частоты, их приведенные массы и интенсивности - характеристики IR спектра.

       1                           1           2           3           4           5
       2        FREQUENCY:       104.70 I     30.68       13.43        9.99        7.70
       3     REDUCED MASS:      1.02197     2.20260     4.20855     4.04406     4.19019
       4     IR INTENSITY:      0.17437     0.00412     0.00002     0.00010     0.00002
    

    Буква "I" после частоты 104.70 означает мнимую единицу. Эта колебательная мода - мнимая. В таком случае говорят о неистинно стационарном состоянии. Требуется более глубокая минимизация. В этом случае переходите к шагу 3, если все в порядке - к шагу 4.

  4. Для более глубокой минимизации установите в групе $statpt параметры "opttol=1e-6" порог, до которого спускаться, и "nstep=200" - максимальное количество шагов спуска.

  5. Откройте програмой GABEDIT ваш гессиановый OUTPUT-файл, визуализируйте самую интенсивную нормальную моду, постройте ИК-спектр.