Сколько отрицательных частот в Hessian?

[kbob]

Нашел, как я предполагаю, переходное состояние в реакции изомеризации кубана.
[cut]16

H 1.408087 1.461905 1.324601
C 0.788578 0.749590 0.780055
H 1.408087 1.461905 -1.324601
C 0.788578 0.749590 -0.780055
H 1.455529 -1.455929 1.483571
C 0.722892 -0.723091 1.179659
H 1.455529 -1.455929 -1.483571
C 0.722892 -0.723091 -1.179659
H -1.389375 1.389759 1.427203
C -0.775771 0.775985 0.771460
H -1.389375 1.389759 -1.427203
C -0.775771 0.775985 -0.771460
H -1.462293 -1.407684 1.324600
C -0.749807 -0.788371 0.780055
H -1.462293 -1.407684 -1.324600
C -0.749807 -0.788371 -0.780055[/cut]

Считаю Hessian в базисе def2-SVP, метод CASSCF(2,2)
Orca дает одну отрицательную частоту

Код: Выделить всё

-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

   0:         0.00 cm**-1
   1:         0.00 cm**-1
   2:         0.00 cm**-1
   3:         0.00 cm**-1
   4:         0.00 cm**-1
   5:         0.00 cm**-1
   6:      -697.33 cm**-1 ***imaginary mode***

NwChem - две

Код: Выделить всё

 ----------------------------------------------------------------------------
 Normal Eigenvalue ||           Projected Infra Red Intensities
  Mode   [cm**-1]  || [atomic units] [(debye/angs)**2] [(KM/mol)] [arbitrary]
 ------ ---------- || -------------- ----------------- ---------- -----------
    1     -708.474 ||    0.000951           0.022         0.927       1.152
    2     -258.583 ||    0.085925           1.982        83.764     104.133

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

[kbob]

Gamess, похоже тоже две отрицательных частоты нашел.

Код: Выделить всё

----- START OF NORMAL MODES FOR -MOLPLT- PROGRAM -----
MODE    1   FREQUENCY= 786.55965 (CM**-1)
MODE    2   FREQUENCY= 310.19132 (CM**-1)
MODE    3   FREQUENCY=  40.31395 (CM**-1)
MODE    4   FREQUENCY=  34.28835 (CM**-1)
MODE    5   FREQUENCY=  32.86577 (CM**-1)
MODE    6   FREQUENCY=   0.00563 (CM**-1)
MODE    7   FREQUENCY=   0.00829 (CM**-1)
MODE    8   FREQUENCY=   0.03885 (CM**-1)
MODE    9   FREQUENCY= 186.42145 (CM**-1)
MODE   10   FREQUENCY= 482.83687 (CM**-1)
MODE   11   FREQUENCY= 601.65688 (CM**-1)
MODE   12   FREQUENCY= 704.70836 (CM**-1)
MODE   13   FREQUENCY= 736.82990 (CM**-1)
MODE   14   FREQUENCY= 792.78120 (CM**-1)
MODE   15   FREQUENCY= 866.89338 (CM**-1)
MODE   16   FREQUENCY= 892.18982 (CM**-1)
MODE   17   FREQUENCY= 893.46888 (CM**-1)
MODE   18   FREQUENCY= 897.35773 (CM**-1)
MODE   19   FREQUENCY= 929.95451 (CM**-1)
MODE   20   FREQUENCY= 931.65597 (CM**-1)
MODE   21   FREQUENCY= 965.36820 (CM**-1)
MODE   22   FREQUENCY=1011.60419 (CM**-1)
MODE   23   FREQUENCY=1018.18109 (CM**-1)
MODE   24   FREQUENCY=1019.11852 (CM**-1)
MODE   25   FREQUENCY=1088.20323 (CM**-1)
MODE   26   FREQUENCY=1115.63153 (CM**-1)
MODE   27   FREQUENCY=1175.97356 (CM**-1)
MODE   28   FREQUENCY=1206.47210 (CM**-1)
MODE   29   FREQUENCY=1257.90377 (CM**-1)
MODE   30   FREQUENCY=1263.55789 (CM**-1)
MODE   31   FREQUENCY=1278.02312 (CM**-1)
MODE   32   FREQUENCY=1279.27487 (CM**-1)
MODE   33   FREQUENCY=1294.64610 (CM**-1)
MODE   34   FREQUENCY=1315.74270 (CM**-1)
MODE   35   FREQUENCY=1355.41931 (CM**-1)
MODE   36   FREQUENCY=1374.48894 (CM**-1)
MODE   37   FREQUENCY=1376.93380 (CM**-1)
MODE   38   FREQUENCY=1387.07794 (CM**-1)
MODE   39   FREQUENCY=1401.59746 (CM**-1)
MODE   40   FREQUENCY=1408.04787 (CM**-1)
MODE   41   FREQUENCY=3229.29606 (CM**-1)
MODE   42   FREQUENCY=3232.14493 (CM**-1)
MODE   43   FREQUENCY=3243.75678 (CM**-1)
MODE   44   FREQUENCY=3244.04059 (CM**-1)
MODE   45   FREQUENCY=3262.77684 (CM**-1)
MODE   46   FREQUENCY=3278.91229 (CM**-1)
MODE   47   FREQUENCY=3371.86016 (CM**-1)
MODE   48   FREQUENCY=3376.09589 (CM**-1)
----- END OF NORMAL MODES FOR -MOLPLT- PROGRAM -----

Незнаю, какой базис в GAMESS соответствует def2-SVP, поэтому посчитано в GBASIS=CCD

[tabaki321]

добрый день. Можно узнать какой метод (HF,DFT,MP2 ...)? А то я увидел только базис. Если скажете метод, смогу попробовать в DALTON2011.

[kbob]

Пардон. Метод CASSCF(2,2), поскольку происходить разрыв связи.
Похоже, что ORCA врет.

[endlesslake]

а энергия и градиент одинаковые ?

[kbob]

а энергия и градиент одинаковые ?

Энергии у меня в ORCA и NwChem всегда отличаются.
Градиенты в этой геометрии близки к нулю - это седловая точка.
В ORCA получается, что это седловая точка 1-го порядка, а в NWCHEM 2-го.

[sanya1024]

Орке доверия мало. Там мноооого всего заявлено, но по-человечески не доделано. Кстати, если Орку еще не обновили до 3-й версии, очень рекомендуется.

[kbob]

ORCA версии 3.0.0. viewtopic.php?f=71&t=101029

У меня к вам 2 вопроса по GAMESS, поскольку у вас большой опыт работы в нем.
1. Как в GAMESS задать базис def2-SVP?
2. Я попробовал расчитать Hessian в базисе def2-TZVP
[cut]$CONTRL SCFTYP=MCSCF RUNTYP=HESSIAN COORD=CART $END
$SYSTEM TIMLIM=600 MWORDS=250 $END
$BASIS GBASIS=KTZVP $END
$GUESS GUESS=HUCKEL $END
$DET NCORE=27 NACT=2 NELS=2 $END
$DATA
Cubane CASSCF(2,2)/TZVP hessian
C1
Hydrogen 1.0 1.408087 1.461905 1.324601
Carbon 6.0 0.788578 0.749590 0.780055
Hydrogen 1.0 1.408087 1.461905 -1.324601
Carbon 6.0 0.788578 0.749590 -0.780055
Hydrogen 1.0 1.455529 -1.455929 1.483571
Carbon 6.0 0.722892 -0.723091 1.179659
Hydrogen 1.0 1.455529 -1.455929 -1.483571
Carbon 6.0 0.722892 -0.723091 -1.179659
Hydrogen 1.0 -1.389375 1.389759 1.427203
Carbon 6.0 -0.775771 0.775985 0.771460
Hydrogen 1.0 -1.389375 1.389759 -1.427203
Carbon 6.0 -0.775771 0.775985 -0.771460
Hydrogen 1.0 -1.462293 -1.407684 1.324600
Carbon 6.0 -0.749807 -0.788371 0.780055
Hydrogen 1.0 -1.462293 -1.407684 -1.324600
Carbon 6.0 -0.749807 -0.788371 -0.780055
$END[/cut]

выдает ошибку

Код: Выделить всё

 ***** ERROR *****
 THIS RUN FOUND    0 CI EIGENVECTORS WITH S= 0.00,
 BUT YOU REQUESTED STATE AVERAGING OF    1 ROOTS.
 PLEASE EXAMINE YOUR CHOICE OF -NSTATE- INPUT DATA.

что-то я не так задал?

[sanya1024]

1. Сторонние базисы задаются через внешний файл. В Гамессе -- нужно в rungms или gms-files найти export EXTBAS= и прописать там путь к файлу с внешним базисом. Плюс ключ extfil=.t. в $BASIS и имя базиса во внешнем файле в ключе GBASIS=. В FireFly по дефолту внешний базис (также включаемый ключом extfil=.t.) живет в файле BASIS.LIB, лежащем рядом с инпутом, или можно в ком. строке после ключа -b задать путь к файлу с внешним базисом. Кстати, для FF есть файл AHLRICHS.LIB (и обновления к нему на сайте Грановского) с def2-базисами.

2. Э... вообще говоря, все неправильно задано. Для MCSCF вообще надо использовать FireFly, в Гамессе Вы результатов не дождетесь, там MCSCF ужасный и неэффективный. Это Вам не Орка -- одно ключевое слово -- и она все сама за Вас сделает (только вопрос, что именно сделает).

В FireFly делается так: сначала надо подготовить вектора стартового приближения. Например, из ХФ. Т.е., сделать ХФ расчет и сохранить группу $VEC из *.dat файла. Потом подготовить и свести CASSCF (скушав в качестве стартовых эти вектора с помощью MOREAD), причем надо отдавать себе отчет в том, состояния какой мультиплетности нас интересуют. У Вас, похоже, бирадикал -- программа сама сгенерила хюккелевские стартовые орбитали и увидела, что с синглетной схемой спаривания ничего не получается, низшее состояние -- триплет. А поскольку Вы понадеялись на дефолты (никогда так не делайте! все проверяйте по мануалу), то у Вас nstate=1 (считать всего одно состояние). Вот программа и ругается: чего, говорит, просишь состояние WITH S= 0.00 (т.е., синглет), да еще усреднять по нему требуешь?

А что Вам Орка-то насчитала? какой мультиплетности состояние поймала она и что за состояние вылезло в других программах? Может, другие программы таки поймали триплет, как оно и должно быть?

Ну, и если нам надо гессиан, то уже сведя CASSCF, надо стартовать расчет гессиана в том же акт. пространстве с теми же настройками, скушав $VEC от сошедшегося CASSCF в качестве начального приближения.

[kbob]

добрый день. Можно узнать какой метод (HF,DFT,MP2 ...)? А то я увидел только базис. Если скажете метод, смогу попробовать в DALTON2011.

Не поделитесть инпутом.

[kbob]

А что Вам Орка-то насчитала? какой мультиплетности состояние поймала она и что за состояние вылезло в других программах? Может, другие программы таки поймали триплет, как оно и должно быть?

В ORCA мультиплетность задается жестко.

Код: Выделить всё

---------------------------------------------
CAS-SCF STATES FOR BLOCK  1 MULT= 1 NROOTS= 1
---------------------------------------------
   
ROOT   0:  E=    -307.1183967793 Eh
      0.70942 [     0]: 20
      0.29058 [     2]: 02

триплетное состояние довольно близко

Код: Выделить всё

---------------------------------------------
CAS-SCF STATES FOR BLOCK  1 MULT= 3 NROOTS= 1
---------------------------------------------
   
ROOT   0:  E=    -307.1108738333 Eh
      1.00000 [     0]: 11

[tabaki321]

К сожалению с CASSCF никогда не сталкивался, поэтому испытываю некоторые затруднения в составлении input-a. Единственное могу подтвердить, что orca порой выдает очень странные результаты даже на DFT.

[sanya1024]

Ну конечно, Орка насильственно заселила орбитали в конфигурации 20, тогда как по жизни тут, видимо, должен быть все-таки бирадикал (конфигурация 11). Следите за руками (в смысле, за программными дефолтами). Разница между бирадикалами (триплетом и синглетом с конфигурацией 11) и замкнутой оболочкой маленькая, а CASSCF оптимизирует орбитали под выбранное состояние, и если состояние выбрано неправильно, будут перекосы. Поэтому очень рекомендуется выбирать сразу несколько состояний (в Вашем случае все еще проще, 2 электрона и 2 орбитали -- надо выбирать все состояния, к-рые получаются на этих орбиталях) и делать по ним усреднение матрицы плотности при CASSCF итерациях. Не уверена, что это можно в Орке, а в FireFly всегда только так и делают (но дефолт -- без усреднения, поэтому внимательно следим).

[kbob]

Ну конечно, Орка насильственно заселила орбитали в конфигурации 20, тогда как по жизни тут, видимо, должен быть все-таки бирадикал (конфигурация 11). Следите за руками (в смысле, за программными дефолтами). Разница между бирадикалами (триплетом и синглетом с конфигурацией 11) и замкнутой оболочкой маленькая, а CASSCF оптимизирует орбитали под выбранное состояние, и если состояние выбрано неправильно, будут перекосы. Поэтому очень рекомендуется выбирать сразу несколько состояний (в Вашем случае все еще проще, 2 электрона и 2 орбитали -- надо выбирать все состояния, к-рые получаются на этих орбиталях) и делать по ним усреднение матрицы плотности при CASSCF итерациях. Не уверена, что это можно в Орке, а в FireFly всегда только так и делают (но дефолт -- без усреднения, поэтому внимательно следим).

ORCA умеет делать усреднение, но это нужно задать руками.
Но тут возникает вопрос: с какми весами усреднять? Вблизи TS возможно с равными весами, а вдоль пути реакции? А в состояниях реагентов и продуктов?
Хотя ваше замечание навело меня на мысль, а не попробовать ли CASSCF(4,4)?

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

[endlesslake]

Ну конечно, Орка насильственно заселила орбитали в конфигурации 20, тогда как по жизни тут, видимо, должен быть все-таки бирадикал (конфигурация 11). Следите за руками (в смысле, за программными дефолтами).

А это разве зависит от программной реализации? Ну кроме симметрии по умолчанию? Где бы про это почитать?

правда для этой системы на 114 ккал выше получается синглетный бирадикал вроде как.

Код: Выделить всё

 STATE   1  ENERGY=      -307.0850259849  S=  0.00  SZ=  0.00  SPACE SYM=A

ALPH|BETA| COEFFICIENT
----|----|------------
 10 | 10 |   0.8426548
 01 | 01 |  -0.5384542

 STATE   2  ENERGY=      -307.0843018456  S=  1.00  SZ=  0.00  SPACE SYM=A

ALPH|BETA| COEFFICIENT
----|----|------------
 10 | 01 |   0.7071068
 01 | 10 |  -0.7071068

 STATE   3  ENERGY=      -306.9038458493  S=  0.00  SZ=  0.00  SPACE SYM=A

ALPH|BETA| COEFFICIENT
----|----|------------
 10 | 01 |   0.7071068
 01 | 10 |   0.7071068

 STATE   4  ENERGY=      -306.8861680845  S=  0.00  SZ=  0.00  SPACE SYM=A

ALPH|BETA| COEFFICIENT
----|----|------------
 01 | 01 |   0.8426548
 10 | 10 |   0.5384542

[endlesslake]

а у орки просто дефолтное смещение для численного дифференцирования неудачное - если задать 0.001 тоже 2 мнимых частоты получаются

Код: Выделить всё

   6:      -762.46 cm**-1 ***imaginary mode***
   7:      -265.88 cm**-1 ***imaginary mode***
   8:       201.69 cm**-1

что не плохо совпадает с аналитическим гессианом в gamess US

Код: Выделить всё

    1     789.092    A1       2.998618    0.026208
     2     288.629    B2       1.163605    1.889622
     3       7.682    B1       4.423560    0.000044
     4       3.216    A2       3.848908    0.000000
     5       1.545    B2       3.820267    0.000071
     6       0.009    A1       6.503901    0.000000
     7       0.017    B2       6.503820    0.000000
     8       7.253    B1       4.560067    0.000040
     9     201.286    A1       1.338664    0.003959

и то и то для одного состояния 20+02.

[kbob]

а у орки просто дефолтное смещение для численного дифференцирования неудачное - если задать 0.0001 тоже 2 мнимых частоты получаются

Код: Выделить всё

   6:      -762.46 cm**-1 ***imaginary mode***
   7:      -265.88 cm**-1 ***imaginary mode***
   8:       201.69 cm**-1

что не плохо совпадает с аналитическим гессианом в gamess US
и то и то для одного состояния 20+02.

Это неплохо совпадает с расчетом в Nwchem (1-й пост).

[sanya1024]

Усреднять надо всегда с равными весами, а жизнь покажет, какое состояние провалится ниже.
Я прошу прощения, слегка народ запутала. Конечно, Гамесс тоже сам по себе бирадикал может не поймать. Дефолтная мультиплетность mult=1, по дефолту pures=.t. (т.е., даже в детерминантных расчетах, для к-рых задается группа $DET, строит чистые по спину состояния для заданной мультиплетности, в данном случае, синглеты). Т.е., если не задать никаких ключей, то Гамесс построит синглетные конфигурации 20, 11 и 02, и совершенно не факт, что 11 окажется ниже, чем 20. Нам же нужны не только синглеты, но и триплеты -- триплетный-то бирадикал (конфигурация 11) уж точно провалится ниже синглетного 11, и весьма вероятно -- что провалится ниже синглета 20 (тут выше показали, что не провалился, чуть-чуть не хватило, но возможно, существует где-то рядом геометрия, где он все-таки провалится).
На пути реакции будет меняться отн. энергия состояний -- главное, не забывать усреднять по ним по всем с равными весами. Если менять веса в усреднении на пути реакции -- это в каком-то смысле будет похоже на смену функционала в DFT. Вы же не станете, скажем, часть точек на ППЭ считать в B3LYP, другую в BHHLYP, а третью в BLYP? вот и тут так же.
Я думаю, переходить к CAS(4,4) рановато. Тем более что в Орке это может быть уже тяжеловато (кол-во конфигураций в акт. пространстве растет с чудовищной скоростью).
И я так и не поняла: "Гамесс" -- это GAMESS-US или FireFly? просто для FireFly есть много полезных ключей, к-рые позволят сильно облегчить жизнь и делать тяжелые расчеты, а в GAMESS-US их нет, и ничего сложного Вы не посчитаете.

[kbob]

а у орки просто дефолтное смещение для численного дифференцирования неудачное - если задать 0.0001 тоже 2 мнимых частоты получаются

Дефолтное значение уменьшено в 50 раз - круто.

Код: Выделить всё

%numgrad
CentralDiff true # (Default) You should use two-sided numerical
# differentiation, but it is possible to switch
# to one-sided numerical differentiation.
DX 0.005
 # Increment in Bohr for the differentiation.
end

[endlesslake]

я опечатался - 0.001 должно быть, и это кажется дефолтное значение для nwchem'a