Релятивистские базисы и порядок связи

[eFo]

Коллеги, научите, пожалуйста, как извлекать из природовского расчета print=+density какую-нибудь информацию, кроме координат критических точек (порядки связей, заряды, etc.)

А нельзя просто добавить +charges+bonds при task=energy?
Что-то типо print=+charges+bonds+density. Тогда в аут-файле всё будет.

[EvgeniX]

Это всё по методу Хиршфельда и Малликена.

[amge]

Я тоже пришел к такому выводу. print=+charges печатает заряды по Хиршфельду, print=+bonds - заряды и порядки связей по Малликену, причем от того, есть там +density или нет, результат, вроде бы, не зависит.

[EvgeniX]

Коллеги, научите, пожалуйста, как извлекать из природовского расчета print=+density какую-нибудь информацию, кроме координат критических точек (порядки связей, заряды, etc.)

QCC Front-End по щелчку на критические точки показывает какие-то параметры, но они такие же, как в выходном файле (вроде a=...., f=...). Тоже был бы рад, если бы кто-нибудь научил

[Лечащий Врач]

Я процентов на 99 уверен, что природа бейдеровские заряды просто не считает. Это очень непросто сделать правильно - имхо только в AIMAll это реализовано как надо, а AIMAll пишет Тодд Кейт, код не раскрывает, а второго Кейта нет

Существует волшебная утилита Molden2AIM, которая якобы делает WFN из молденовских файлов. Если кто-нибудь попробует это проделать с выводом Природы и положит сюда результирующий WFN (скажем, для циклопропана) буду очень признателен, потому что сам природой считать не умею, но если окажется, что все-таки делать WFN из нее можно, стану осваивать

[EvgeniX]

Расчёт циклопропана.

Файлы:
Cyclopropane_PBE_L2_o.in - Входной файл оптимизации структуры циклопропана методом PBE/L2;
Cyclopropane_PBE_L2_o.out - Выходной файл оптимизации;
Cyclopropane_PBE_L2_h.in - Входной файл решения колебательной задачи;
Cyclopropane_PBE_L2_h.out - Выходной файл решения колебательной задачи;
Cyclopropane_PBE_L2_b.in - Входной файл для проведения бейдеровского анализа средствами ПРИРОДы;
Cyclopropane_PBE_L2_b.out - Выходной файл бейдеровского анализа.

Секция выходного файла ПРИРОДы, посвящнная бейдеровскому анализу:

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

 Critical Points of the Density
 nx=43, ny=42, nz=37
...........................................................................
 Np=15
 Np=10

 r=  0.989414, -0.752189, -0.579218 f=0.278034 h=-0.73757,-0.71773, 0.58212
 r= -1.146119, -0.480746, -0.579208 f=0.278033 h=-0.73757,-0.71773, 0.58212
 r=  0.156705,  1.232932, -0.579207 f=0.278033 h=-0.73756,-0.71773, 0.58212
 r= -0.058852, -0.463158, -0.000006 f=0.232729 h=-0.43872,-0.28775, 0.41556
 r=  0.000016, -0.000016, -0.000003 f=0.193488 h=-0.34329, 0.27008, 0.27011
 r=  0.430522,  0.180593, -0.000004 f=0.232726 h=-0.43871,-0.28774, 0.41556
 r= -0.371686,  0.282548,  0.000001 f=0.232741 h=-0.43875,-0.28779, 0.41558
 r=  0.989442, -0.752123,  0.579213 f=0.278034 h=-0.73757,-0.71773, 0.58212
 r= -1.146085, -0.480821,  0.579208 f=0.278033 h=-0.73757,-0.71773, 0.58212
 r=  0.156642,  1.232942,  0.579210 f=0.278033 h=-0.73757,-0.71773, 0.58212

/* Координаты атомов */
 Atomic Coordinates:
   6  -0.80306417  -0.33688018  -0.00000080
   6   0.10976859   0.86390871   0.00000195
   6   0.69328701  -0.52704722  -0.00000788
   1  -1.34857118  -0.56579393   0.91475564
   1  -1.34862971  -0.56566679  -0.91475409
   1   0.18430768   1.45078647   0.91475683
   1   0.18441717   1.45076883  -0.91475570
   1   1.16421771  -0.88509446  -0.91476542
   1   1.16426694  -0.88498139   0.91476953

/* Координаты критических точек */
   0   0.98941425  -0.75218894  -0.57921759
   0  -1.14611935  -0.48074608  -0.57920783
   0   0.15670492   1.23293192  -0.57920743
   0  -0.05885195  -0.46315801  -0.00000560
   0   0.00001564  -0.00001561  -0.00000275
   0   0.43052246   0.18059321  -0.00000374
   0  -0.37168553   0.28254813   0.00000136
   0   0.98944206  -0.75212306   0.57921338
   0  -1.14608477  -0.48082066   0.57920817
   0   0.15664208   1.23294194   0.57920988
 #

/* Время, затраченное на анализ электронной плотности (по Бейдеру) */
 time: Density                          :    29.3  97.8%    30.0

[EvgeniX]

Расчёт циклопропана, продолжение.

Файлы:
Cyclopropane_PBE_L2_wf.in - Входной файл для получения волновой функции в формате Molden;
Cyclopropane_PBE_L2_wf.out - Выходной файл c волновой функции в формате Molden;
Cyclopropane_PBE_L2_wf.molden - Файл Molden с волновой функцией, полученный из Cyclopropane_PBE_L2_wf.out;
Cyclopropane_PBE_L2_wf.wfn - WFN-файл, полученный из Cyclopropane_PBE_L2_wf.molden с помощью Molden2AIM.

В программе Molden2AIM указано, что в настоящее время тестируется на совместимость с программами StoBe, Columbus и PRIRODA!

[amge]

Существует волшебная утилита Molden2AIM, которая якобы делает WFN из молденовских файлов. Если кто-нибудь попробует это проделать с выводом Природы и положит сюда результирующий WFN (скажем, для циклопропана) буду очень признателен, потому что сам природой считать не умею, но если окажется, что все-таки делать WFN из нее можно, стану осваивать

Я попробовл эту утилиту: действительно, она делает WFN и, вроде бы, правильно. Вот bash-скрипт, облегчающий ее использование для Природы:

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

#!/bin/sh

#set -x
if [ ! -n "$1" ]; then
  echo "Usage: `basename $0` Priroda_output_filename"
  echo "(print=+vectors+molden must be in Priroda input)"
  exit
fi

m2a=/usr/local/bin/m2a
if [ ! -x $m2a ]; then
  echo No Molden2AIM executable
  exit
fi

for file in $*; do
  TEMP=`tempfile -d .`
  #perl -ne 's/mos>// && print' $file > $TEMP
  grep 'mos>' $file | sed -e 's/mos>//' > $TEMP
  if [ ! -z $TEMP ]; then
    echo -ne "\n$TEMP\n\n\n" | $m2a > /dev/null
  fi
  rm $TEMP
  if [ -f $TEMP.wfn ]; then
    name=`basename $file .out`
    mv $TEMP.wfn $name.wfn
  fi
done

Upd. О! Пока я писал скрипт, EvgeniX меня опередил

[Лечащий Врач]

Теперь моя очередь

Я быстренько посчитал циклопропан в Гауссиане PBEPBE/cc-pVTZ (настолько близко к природовскому PBE/L2, насколько возможно).

В приложении WFN, sum-файлы программы AIMAll для обоих расчетов, а также сводка КТ в xls (мой собственный скрипт генерирует, так читать данные по КТ сильно удобнее). Обратите внимание, что плотность и лапласиан там приведены в e*A^-3! Энергии в конце - по корреляции Эспинозы/Лекомта.

Краткий анализ:

• параметры КТ внешне идентичны (отличие менее 1%, даже для лапласиана точек (3,-1)) — уже неплохо
• ОДНАКО параметры КТ в AIMAll для WFN из природы отличаются от рассчитанных собственным анализом природы: для связи С-С природа насчитывает rho(r)=0.2327, а AIMAll по сгенеренному WFN — 0.2352. Собственные значения гессиана тоже другие. Откуда разница?
• с интегральными характеристиками атомов сложнее — интегрирование WFN из природы дает не вполне корректный результат для зарядов атомов углерода. Согласно фундаментальной работе Бейдера по углеводородам, у углерода в циклопропане должен быть небольшой отрицательный заряд. Гауссиан дает -0.055 е, а в природе получается +0.008 — различие большое.
• Делокализационные индексы (табличка в конце, величина DI(A,B)) врут на сопоставимую с зарядами величину: в Гауссиане для связи C-C 1.0005 (корректное значение), а в природе 0.9944
• для корректного анализа энергий атомов нужно знать итоговую энергию и вириал (V/T). Их нужно вписывать в WFN руками (я не стал). Но ошибка в 0.5% для энергии атома - это очень много, а относительные ошибки энергии обычно имеют тот же порядок, что и ошибки зарядов.

Похоже, что природовский алгоритм поиска КТ совсем дурной. Если бы ошибки в КТ отвечали ошибкам округления, ошибки интегрирования атомов (лапласианы) имхо были бы больше.

Но делокализационные индексы, если иметь это в виду и сравнивать только природовские расчеты с природовскими же, применять, похоже, возможно.

С зарядами жалко — для моих задач ошибка в 0.06е очень велика. Например, сопоставимое изменение возникает за счет сильного аномерного эффекта. Правда, возможно, что наблюдаемые различия в интегральных характеристиках вызваны реализацией PBE в природе. Интересно, что будет в riMP2?

[EvgeniX]

Очень интересно

[*] ОДНАКО параметры КТ в AIMAll отличаются от тех, которые считает природа: для связи С-С природа насчитывает rho(r)=0.2327, а AIMAll по сгенеренному WFN — 0.2352. Собственные значения гессиана тоже другие. Откуда разница?

Это сравнивается Gaussian(WFN)+AIMAll и непосредственно анализ в ПРИРОДе (параметр "f=...")?

[*]с интегральными характеристиками атомов сложнее — природа выдает не вполне корректный результат для зарядов атомов углерода. Согласно фундаментальной работе Бейдера по углеводородам, у углерода в циклопропане должен быть небольшой отрицательный заряд. Гауссиан дает -0.055 е, а в природе получается +0.008 — различие большое.

Это ведь уже ПРИРОДА(wfn)+AIMAll?

RI-MP2/L2-расчёты:

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

 Critical Points of the Density
 nx=41, ny=43, nz=37
...........................................................................
 Np=10
 Np=10
 Np=10
 r=  0.980373, -0.745329, -0.577464 f=0.292137 h=-0.80948,-0.79224, 0.45698
 r= -1.135681, -0.476332, -0.577444 f=0.292134 h=-0.80946,-0.79223, 0.45699
 r=  0.155301,  1.221661, -0.577448 f=0.292134 h=-0.80946,-0.79223, 0.45699
 r= -0.058946, -0.463852, -0.000011 f=0.241511 h=-0.47926,-0.31399, 0.34022
 r=  0.431136,  0.180843, -0.000007 f=0.241499 h=-0.47923,-0.31395, 0.34022
 r=  0.000021, -0.000024, -0.000004 f=0.200422 h=-0.37972, 0.25893, 0.25897
 r= -0.372246,  0.282959,  0.000006 f=0.241526 h=-0.47930,-0.31404, 0.34022
 r=  0.980404, -0.745220,  0.577464 f=0.292137 h=-0.80947,-0.79223, 0.45699
 r= -1.135617, -0.476480,  0.577444 f=0.292134 h=-0.80946,-0.79223, 0.45699
 r=  0.155213,  1.221686,  0.577445 f=0.292135 h=-0.80946,-0.79223, 0.45699
 Atomic Coordinates:
   6  -0.79915673  -0.33524314  -0.00000038
   6   0.10921476   0.85969335   0.00000321
   6   0.68989819  -0.52449787  -0.00001290
   1  -1.33393598  -0.55973864   0.90948498
   1  -1.33404446  -0.55948681  -0.90948433
   1   0.18232210   1.43506062   0.90948445
   1   0.18247756   1.43501782  -0.90949380
   1   1.15158284  -0.87549597  -0.90951067
   1   1.15164173  -0.87530941   0.90952945
   0   0.98037250  -0.74532904  -0.57746420
   0  -1.13568076  -0.47633163  -0.57744402
   0   0.15530089   1.22166098  -0.57744755
   0  -0.05894593  -0.46385172  -0.00001104
   0   0.43113633   0.18084286  -0.00000660
   0   0.00002078  -0.00002415  -0.00000434
   0  -0.37224584   0.28295907   0.00000552
   0   0.98040380  -0.74522021   0.57746430
   0  -1.13561661  -0.47648007   0.57744423
   0   0.15521331   1.22168637   0.57744486
 #
 time: Density                          :    95.5  [    47.7    47.7 ] 100.0%    47.7

upd. WFN перезалил

[Лечащий Врач]

Очень интересно
Это сравнивается Gaussian(WFN)+AIMAll и непосредственно анализ в ПРИРОДе (параметр "f=...")?

Это ведь уже ПРИРОДА(wfn)+AIMAll?

Да, все правильно, поправил исходный постинг для большей однозначности. А WFN, который вы прикрепили - это не WFN

[EvgeniX]

Перезалил.

[Лечащий Врач]

КТ поменялись, заряды атомов улерода стали +0.09

Видимо, природовские density-fitted базисы неидеальны для изучения электронной плотности Или не могут быть использованы напрямую. Но тут уже идет теория, которой я не знаю.

Скорее всего, впрочем, для межмол взаимодействия и оценки переноса заряда все будет работать нормально (опыт расчета в плоских волнах и мультипольного уточнения подсказывает, что если ковалентные связи описаны неидеально, то межмол все равно может быть ОК).

[EvgeniX]

Ошибку в Molden2AIM нашли

строки1746-1757: c MaxL

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

if(ng.gt.0)then
MaxL=4
else if(nf.gt.0)then
MaxL=3
else if(nf.gt.0)then
MaxL=2
else if(nf.gt.0)then
MaxL=1
else
MaxL=0
end if

вместо nf'ов должно быть что-то в порядке ng.gt.0, nf.gt.0, nd.gt.0, np.gt.0
Как это повлияет на работу программы пока не понятно Используется в секции проверки ВФ (chkbstyp).
Может автору написать?

[amge]

Ошибку в Molden2AIM нашли Может автору написать?

Написать, конечно. Заодно можно сообщить, что с Природой работает.

Кстати, в связи с тем, что автор Molden2AIM - китаец, вспомнил, что меня одно время донимали какие-то китайцы (а может, это один и тот же был) с просьбой поделиться Природой. Приходилось посылать (мейл Лайкова). Теперь можно давать ссылку на SI: Природа оттуда хоть и непараллельная, но вполне полнофункциональная, только, чтобы ею обычные расчеты делать, прежние базисы нужны, идущий в комплекте basis.in для этого не годится.

[EvgeniX]

Написать, конечно. Заодно можно сообщить, что с Природой работает.

С автором списался, это действительно опечатка Обещал исправить в ближайшее время.

[EvgeniX]

В Molden2AIM обнаружили ещё пару мест (более важных, где происходит расчёт), где реализованные формулы не совпадают с формулами, указанными в комментариях (в одном месте степень не та и т.п.). Пока знаний не хватает разобраться, то ли формулы не верны, то ли комментарии к ним.

[amge]

Может, стоит опять написать автору, указав на эти несоответствия. Пусть разбирается.

[surius]

да здравствует некро... ну в общем:

ОДНАКО параметры КТ в AIMAll для WFN из природы отличаются от рассчитанных собственным анализом природы: для связи С-С природа насчитывает rho(r)=0.2327, а AIMAll по сгенеренному WFN — 0.2352. Собственные значения гессиана тоже другие. Откуда разница?

Некоторое время назад натыкался на статейку, в которой всторы показали что практически все выходные данные Бадеровского анализа сильно 
зависят от базисного набора. Особенно ярко такое наблюдается для ковалентных связей. А резюме было такое, что зачастую даже 
6-311++g(d,p) не достаточно, а авторы рекомендуют пользовать сс-pvtz. Так что может быть эта разница набежала за счет различия в базисных наборах?!

А в качестве ремарки, я тут сравнивал поведение различных зарядов  и нашел такую штука - наименее чувстветельным к базисному набору является анализ по 
Хиршфельду. Как только базис перерос 6-31+g, то они (заряды) стоят как вкопанные. А еще иногда удается установить прямую корреляцию между анализом 
зарядов по Бадеру и Хиршфельду. В вертикальных процессах они видут себя одинаково.

[Лечащий Врач]

Некоторое время назад натыкался на статейку, в которой всторы показали что практически все выходные данные Бадеровского анализа сильно 
зависят от базисного набора. Особенно ярко такое наблюдается для ковалентных связей. А резюме было такое, что зачастую даже 
6-311++g(d,p) не достаточно, а авторы рекомендуют пользовать сс-pvtz. Так что может быть эта разница набежала за счет различия в базисных наборах?!

Она очевидно набежала за счет разницы в базисах Тем не менее, L2 в природе по размеру не уступает cc-pvtz, и вообще на него очень похож (и в гауссиане я считал как раз в последнем). Но мы не знаем, что именнно она выдает в molden и в каком виде туда попадает фирменный природовский дополнительный базис. Если ни в каком, то чего ждать?

6-311++g(d,p) достаточно почти всегда, но иногда Попловские базисы на коротких связях между элементами третьего периода приводят к ложным неядерным аттракторам. Но если у вас нет связей P=S или им подобных, то можно этого не опасаться.

Вообще параметры критических сточек слабо зависят от базиса, заряды - заметнее. Еще они зависят от того, каким методом считать.

А в качестве ремарки, я тут сравнивал поведение различных зарядов  и нашел такую штука - наименее чувстветельным к базисному набору является анализ по 
Хиршфельду. Как только базис перерос 6-31+g, то они (заряды) стоят как вкопанные. А еще иногда удается установить прямую корреляцию между анализом 
зарядов по Бадеру и Хиршфельду. В вертикальных процессах они видут себя одинаково.

Не вы первый Особенно это хорошо работает в координационных соединениях. Вне работы просто не найду статьи.