Не сходится Хартри-Фок, кто виноват и что делать?

[leonidas]

Коллеги,

может кто сталкивался с подобным. Считаю энергию взаимодействия для комплекса бензол-ксенон, hf+ccsd(t), без ECP.
Считаю в Molpro, появляется ошибка
?Error: RHF not converged. This error exit can be avoided using the IGNORE_ERROR option on the ORBITAL directive

Тут рецепты на эту болезнь
https://www.molpro.net/info/2015.1/doc/ ... de177.html

Перепробовал все варианты кроме последнего, использовать MCSCF, с которым у меня нет опыта.
Если же сделать как рекоммендует аутпут Molpro, т.е. IGNORE_ERROR, тогда выходит хрень на уровне CCSD
UNREASONABLE NORM. CALCULATION STOPPED

Еще варианты попробовать другой пакет, например Gaussian, имеет ли смысл?
Уже опускаются руки, вся надежда на Вас...

[chemigor]

Я в этом пакете не считал, да и в гауссиане в прошлом тысячелетии считал, но сам программировал Хартри-Фок и занимался его сходимостью, ибо диссер у меня был об этом: http://www.dissercat.com/content/algebr ... voi-khimii

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

ИМХО, если можете чуть-чуть подправить свой базис, то сходимость должна улучшится.

[endlesslake]

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

С molpro дело не имел, не может быть такого, что вы на ксеноне задали базис без остовных функций, а ECP не задали ? в Gamess'е это была распространенная причина подобных проблем.
не соответствие выбраного учета скалярно-релятивистских тому для чего был оптимизирован базис, по идее, тоже может проблемы создавать.

[Гесс]

Расчет как HF падает аналогично?
Расчет с меньшим базисом (каким?) сводится но чтение его в большем базисе не помогает?
Так а как шаги сведения выглядят, оно осцилирует? Или ему шагов не хватило?
DF не используется?
Какие еще конвержеры доступны в Молпро кроме ДИИСа? (я так понимаю у вас он). Если есть варианты я бы начал с этого пункта. Ньютон-Рапсон, КДИИС, СОСЦФ?

[leonidas]

Большое спасибо всем за ответы!

Забыл уточнить что не сходится только энергия димера, т.е. Eab,
Ea и Eb худо-бедно сходятся. Точнее Eab сходится к неправильному значению
при использовании сдвига уровней (level shift), при этом энергия взаимодействия
Eint=Eab-Ea-Eb выходит абсурдная.

Ниже детали, если кому-то будет интересно.

1. Я писал в соседней ветке, что small-core ECP для комплексов
с Xe врёт - даёт слишком мелкий потенциал.
Это сподвигло меня на то, чтобы посчитать явным способом корреляцию остовных электронов.
Да, долго и дорого, но считает. Здесь проблем со сходимостью нет в базисах типа tight, acvXz.
В данном случае речь о ECP не идёт вообще!

2. Если не использую ECP, тогда нужно каким-то другим способом учитывать релятивисткие поправки.
Для подобных комплексов около 90% этих поправок даёт
Scalar-relativistic Douglas-Kroll-Hess (DKH) Hamiltonian.

Причем для подобных вычислений нужно использовать базисы типа decontracted.
Вот с ними у меня и возникают проблемы со сходимостью для Eab!

3. Базисы (не ECP!) которые я использую для ксенона достаточно новые и пока что мало где использовались.
Мое мнение что они немного "самопальные", но у меня нет времени создавать и тестировать свои.
Если нужно могу дать ссылки. Их тестировали для Xe-Xe, с которым проблем не возникало.

4. Если считать HF отдельно, тогда считается до конца, при этом остается Warning в аутпуте. Но непонятно какой смысл
этих вычислений. DKH нужно считать именно с CCSD(T).
Вы это имели ввиду?

5. Если не использовать level shifts, тогда энергии осциллируют, в аттачменте аутпут.

6. В Molpro доступны только 2 конвержера, DIIS и его вариация KAIN (Krylov-subspace accelerated inexact Newton).
В каких пакетах доступны упомянутые Вами? Может имеет смысл их попробовать?

7. Начальное решение с меньшим базисом cc-pVDZ для C и H, adzp для Xe не помогает для Eab.

8. Существенная особенность моих вычислений - использование дополнительных функции типа midbond, которые улучшают сходимость
Eint когда используются маленькие орбитальные базисы типа avdz. Без этих функций midbond всё сходится, иногда с помощью начального вычисления с меньшим базисом, но при этом результаты не надежные. Нужно или брать базисы типа av5z, что почти невозможно для моего комплекса, или брать базисы avdz+midbond.

CBS(X,X+1) интерполяция? Вещь тоже ненадежная для X=2, но это отдельная тема.

9. Геометрия Re, т.е. координаты глобального минимума ППЭ (более-менее).

10. DF не использую, так как такой возможности нет в Molpro. Имеет ли смысл попробоать DF-CCSD(T) с DKH в Орке?
А может даже DLPNO с DKH? При этом придется брать левые auxiliary базисы, так как для Xe ещё нет базисов типа mp2fit и jkfit,
что тоже даст дополнительную ошибку.

Такие дела...

[chemigor]

Причем для подобных вычислений нужно использовать базисы типа decontracted.
Вот с ними у меня и возникают проблемы со сходимостью для Eab!

я не спец в этом пакете, но что-то мне подсказывает, что деконстрейнед - это не отртогональный базис. Если так, то понятно Ваша решалка подправляет (доортогонализует этот базис) перед началом самосогласования и даже перед первой Хартри итерации и вот тут может получаться очень "кривые" функции. Я когда сам такое же писал, то обычно сами базисные функции как-то визуализировал, и, если они на мой "химический" взгляд казались не разумными, тут же бежал и правил базис. Возможно такой метод и Вам поможет (если он конечно реализуем).

Еще вариант - посмотреть на решение (визуализацию) Хартри решения, если оно кривое, то сходится тоже не обязано.

При адекватном базисе и 10 миллионах базисных функций у меня сходимость самосогласования бывала меньше 10 итераций (примерно до 1е-7 по невязке решения), но я в свое время именно доказывал, чтоб такое происходило.

[endlesslake]

Причем для подобных вычислений нужно использовать базисы типа decontracted.

этот тезис кажется немного сомнительным. То есть понятно, что нужно использовать базисы recontracted для используемого учета релятивизма, а зачем нужен именно decontracted для энергии нековалентного взаимодействия, мне как то не понятно.
Или предполагается, что таким образом как то лучше учтется core-valence и core-core correlation ?

Кстати, я из лога не смог сходу понять, проводит ли molpro чистку базиса от линейно зависимых функций? Как мне кажется, без этого, при таких способах составления базиса (mid-bond и\или decontraction) проблемы со сходимостью будут все время возникать.

[alien308]

Насколько понимаю мольпро имеет свой, довольно гибкий язык для составления заданий. Возможно на нём можно посчитать молекулы по отдельности и использовать их орбитали для стартовой конфигурации "fragment superposition".

[Гесс]

Перечисленные конвержеры есть в орке. Вгауссиане есть ньютонрафсон под именем scf=xc. В обоих пакетах это самый мадленный метод позволяющий сводить самые плохие случаи. Плюс у орки есть советы что еще с ней делать при несводимости. Я слышал мнение что орка сводится медленнее но надежнее чем гауссиан. Впрочем возможно это касалось лишь дфт. Кроме дкх в орке для релятивизма есть зора совместимо ли это с длпно не в курсе. Ауксиларники для нестандартных сетов орка умеет создавать сама и эта процедура была переработана в 4.0 релизе. Если не устраивает дз-тз экстраполяция а 5з слишком дорого - тз-кз? Для улучшения коркор и пр взаимодействий можно вместо традициногл даннинга ccpvxz брать ccpcvxz или его w weighted аналог но оно заметно дороже. Сорри я с телефона.

[YuraM]

Интересно, а зачем использовать какие-то необычные (по-крайней мере для обычных расчетов) базисы? Для углерода есть cc-pwCVnZ, для водорода cc-pVnZ, для ксенона есть вроде cc-pwCVnZ-PP или cc-pwCVnZ-DKH которые как раз созданы для описания корреляции подвалетных оболочек, или хочется коррелировать уж совсем остовные электроны, например 1s 2s2p 3s3p у ксенона? Не, они конечно какой-то вклад дадут, но что-то мне подсказывает - очень маленький, те нужна какая-то супер-точность?

http://www.grant-hill.group.shef.ac.uk/ccrepo/xenon/

[leonidas]

Попробовал посчитать тоже саме в Gaussian и вообще не возникло проблем со сходимостью, кроме одного случая
когда энергию Ea(benzene) надо было считать в 2 шага, т.е. сначала в меньшем базисе. При этом использовал
конвержер по умолчанию, т.е. EDIIS. Кстати Гесс имел ввиду видимо XQC? Я в лёкгом недоразумении, так как всегда был уверен в Molpro.

Отвечаю на комментарии.

1. Использование функции типа midbond является де-факто стандартом для вычисления энергии взаимодействия
малых и средних комплексов супермолекулярными методами и SAPT. Исключение - методы R12/F12 для которых это не нужно.

2. Действительно recontracted и decontracted это не одно и тоже. Просто с первыми возникло больше проблем
со сходимостью в Molpro, поэтому решил использовать вторые. Для Xe есть базисы типа avXz-dk и именно их
буду использовать.

3. Для C и H я использовал стандартные базисы, про PP я писал выше почему я от него отказался.
Кстати, спасибо за ссылку, гораздо более удобно чем EMSL, который постоянно лежит.

4. Насколько я знаю Molpro проверяет линейную зависимость и выдает ошибку в случае чего.

[Гесс]

Да, сорри, Гесс имел в виду scf=qc,
scf=xqc начинает со стандартного diis и если не сходится за половину scf шагов то переходит на ньютонрафсона, scf=qc сразу начинает с нбютонР. xqc - это мой традиционный кейворд, он дает системе шанс сойтись по хорошему и если не получается то переходит к поплохому. Разумеется есть ряд недостатков, например: диису нехватило всего пары шагов когда в бой вступил qc (получается намного дольше чем без этого кейворда) или диис сошел с ума и передал qc намного худшее начальное приближение чем было вначале (врезультате qc не сводится или сводится нетуда и это занимает недели).

Кстати а почему R12/F12 ненужен мидбонд?

[leonidas]

Не будучи эспертом по r12/f12, насколько я понимаю вся эта теория r12/f12 включения парных функций взаимодействия электронов, т.е. explicit correlation, была затеяна для более быстрой сходимости. Как результат, базис avtz для димеров уже дает значения близкие к пределу.

[leonidas]

Тут есть ещё такая особенность, чтобы найти релятивисткую поправку явно (речь о DKH2), нужно посчитать энергии с и без этой поправки
в той же самом базисе, т.е. avxz-dk. Можно конечно вычисления без поправки считать в базисе avxz, что работает для легких атомов (1-2 ряд),
но для более тяжелых уже даёт большую разницу.

Gaussian справился с релятивисткими вычислениями, а вот с нерелятивисткими дал нефизичиские значения.
Пустил с опцией stable и действительно даёт сообщение:
"The wavefunction has an RHF -> UHF instability."

т.е. SCF сходится к седлу или локальному минимуму. Что делать, пробовать другие конвержеры, попробовать guess=mix?

Полезные ссылки по теме:
http://www.ccl.net/cca/documents/dyoung ... verge.html

http://www.cup.uni-muenchen.de/ch/compc ... ility.html

[Гесс]

Пустил с опцией stable и действительно даёт сообщение:
"The wavefunction has an RHF -> UHF instability."

т.е. SCF сходится к седлу или локальному минимуму. Что делать, пробовать другие конвержеры, попробовать guess=mix?

stable=opt
насколько я могу судить - сцф сводится к минимуму, но неправильному минимуму (решению).
предположительно (по хрустальному шару) - вы считаете rhf а оно утверждает что существует более низкое uhf решение. (возможно кстати под этим uhf решением есть и еще более низкий rhf но тут в хрустальном шаре все затянуто кровью.

[YuraM]

May be I am wrong, but in Gaussian if you specify stable=opt and UHF as a calculation method, an attempt to find the stable function will be taken, and the stable wavefunction can be saved (and later read by another calculation). However, if you specify RHF + stable=opt the claculation says there is RHF UHF instability and that is it.

[Гесс]

тут не помню, возможно зависит от версии гауссина, вроде в последней версии stable=opt будет искать по uhf если это потребуется. Незнаю подавляется ли это прямым указание rhf. коль скоро цена ошибки - цена синглпоинта - не вижу проблемы проверить методом тыка.

[leonidas]

У меня опция stable=opt не работает, работает только stable. При этом Gaussian игнорирует ccsd(t) в инпуте, что в принципе в данном случае и не нужно.

Хрустальный шар вдруг осветился ярким светом и выдал файл с аутпутем, а потом погас .

[YuraM]

I am not the best coupled cluster expert, but there are a few things in your output that worry me. First, why option stable=opt would not work in your case? I also used G09 and it worked. My wild guess is that it is available for HF and DFT, Could you try to run HF calculation with this option? Then, please, explicitly specify UHF. I think by default RHF is used for the singlets, which is confirmed in your input: RCCSD(T)-FC. When you get your wave function converged you can read it and use for your ccsd(t)

Then, if I am not mistaken, you wanted to correlate the subvalence electrons, but your CCSD(T) uses frozen core, please check the frozen core defaults with the basis sets you use.

[leonidas]

Опция stable=opt работает с самим HF и вроде делает то же самое, что stable, но при этом даёт более развернутый аутпут.

Конечно, было бы здорово сразу считать всё с остовными электронами, но на практике обычно считают в приближении FC в маскимально возможном базисе, и потом, в меньшем уже поправки типа core-valence и core-core, а также типа T-(T) и (Q)-T для CC.
См. напр. 10.1063/1.4994267