RedCar failed в ходе Relaxed PES

[Phoebe12]

Здравствуйте, уважаемые форумчане!

Я пытаюсь сделать Relaxed PES в Gaussian с целью оценки возможной области переходного состояния и дальнейшего поиска этого TS. Сканирование выполняется для реакции переноса протона с аминокислотного остатка тирозина на молекулу воды при условии, что в окружении донора и акцептора есть ещё несколько молекул воды и аминокислотных остатков. Исходные координаты всех атомов взяты из pdb-файла белка (кроме водородов, которые, как известно, в pdb обычно не присутствуют), при этом для имитации ограничений движения, которые имеют место в белке, координаты некоторых атомов заморожены.

Исходная структура:

Изначально я пыталась сделать сканирование на основе B3LYP/6-311++G, не оптимизируя заранее систему "протонированный донор+акцептор+их окружение", но эта попытка не дала ничего хорошего. На каждой итерации сканирования выполнялось 276 из 276 циклов оптимизации, геометрия при этом не оптимизировалась, и система переходила к следующему шагу сканирования. Затем я попробовала предварительно оптимизировать систему реагентов сначала на основе HF/STO-3G, после этого B3LYP/STO-3G и, в конце концов, B3LYP/6-311++G. Такая последовательная тройная оптимизация дала следующую структуру:

Оптимизированная структура:

С этой геометрии я начала Relaxed PES на основе B3LYP/6-31G. На первом шаге сканирования выполнились 3 из 4 критериев сходимости. На этапе расчёты крашатся, поскольку "RedCar failed".

Прилагаю текст input`a:

%chk=scan.chk
%mem=1Gb
%nproc=8
# B3LYP/6-31g Opt=modredundant

From Tyr33 to H2O proton transfer relaxed scan

1 2
C 6.1010170 0.2495100 -2.1610350
C 5.4804120 0.3468000 -0.7126070
C 3.9988290 -0.0099230 -0.4721280
C 3.5463770 -1.3626340 -0.4060470
C 3.0732370 1.0305510 -0.1455200
C 2.2555000 -1.6688730 0.0334070
C 1.7762890 0.7540220 0.2939100
C 1.3429220 -0.6167110 0.4210560
O 0.1224110 -0.9885200 0.9377710
H 5.6146720 1.3855610 -0.3670450
H 6.0883850 -0.2948770 -0.0477300
H 4.2391960 -2.1765020 -0.6641480
H 3.4109830 2.0754790 -0.1999360
H 1.9177870 -2.7126000 0.1148030
H 1.1061400 1.5762530 0.5872400
H 5.4241940 0.6592200 -2.9246510
C -3.4231570 -2.4289420 -1.8443690
C -2.9908250 -2.3490590 -0.3503180
O -4.1806910 -2.6371660 0.4640180
H -2.5811190 -1.3319250 -0.1357790
H -2.1678070 -3.0785140 -0.1638980
H -3.7831920 -2.6648210 1.4127320
H -4.2759960 -1.7562450 -2.0188850
C -5.1527540 1.7537620 -1.3099380
C -3.7802140 2.3319960 -1.7743290
O -2.7098460 1.7073610 -0.9043060
H -3.7520470 3.4391640 -1.6532500
H -3.5662830 2.0935930 -2.8386880
H -1.8717520 1.9811170 -1.4259710
H -5.1770640 0.6810880 -1.5596990
O -1.4954400 0.7207380 1.4359900
H -1.0242740 1.6424000 1.7438600
H -2.0330750 1.0140650 0.5787300
O -0.2921400 3.1085020 2.1428470
H -0.0078820 3.0602730 3.1219710
H -1.0135140 3.8271200 2.2158900
O 0.3322130 -1.0086940 3.1526120
H 1.1356870 -1.6347700 3.2729730
H -0.4095040 -1.6689810 3.4092160
H 6.9563810 0.8855880 -2.0680000
H 6.3646140 -0.7579940 -2.4067020
H -5.0571680 1.8201150 -0.2462840
H -6.0991910 2.1571470 -1.6039440
H -2.6198000 -2.1363580 -2.4877160
H -3.7303840 -3.4281220 -2.0727020
H -0.6658340 -0.0572290 1.1636750

1 F
17 F
24 F
31 F
34 F
37 F
40 F
41 F
42 F
43 F
44 F
45 F
B 46 31 S 4 -0.0674

Соответствующие строки output`а:

Item Value Threshold Converged?
Maximum Force 0.000004 0.000450 YES
RMS Force 0.000001 0.000300 YES
Maximum Displacement 0.002380 0.001800 NO
RMS Displacement 0.000360 0.001200 YES
Predicted change in Energy=-3.389168D-08
Optimization completed on the basis of negligible forces.
-- Stationary point found.

..... (пропускаю часть output-файла, где приведена таблица Optimized Parameters)

Iteration 1 RMS(Cart)= 0.00476819 RMS(Int)= 0.00581262
New curvilinear step failed, DQL= 6.89D-02 SP=-7.94D-03.
RedCar failed.

Кто-нибудь может подсказать, в чём причина и как это можно исправить/обойти? У меня сложилось впечатление, что проблема кроется либо в самой исходной геометрии, с которой я начинаю сканирование, либо в совокупном неверном задании ограничений и параметра сканирования (если вообще не во всём сразу). Заранее благодарю!

[Гесс]

На мой взгляд схема "HF/STO-3G, после этого B3LYP/STO-3G и, в конце концов, B3LYP/6-311++G" довольно длинна, а 6-311++G с одной стороны слишком велик, (расчет идет весьма медленно) а с другой - не содержит поляризационных функций. Я бы предложил вместо всей этой схемы поробовать чтото типа B3LYP/6-31+G* и переходить к большему базису лишь когда все будет выяснено в этом. Проблемы с curvilinear step часто связаны с внутренними координатами, возможно в подробном аутпуте (при указании в инпуте P - я имею в виду "#P B3LYP...") - будет видно какие координаты проблемные. Кардинальным решением проблемы может быть использование декартовых координат для оптимизации (Opt=Cartesian,Modredundant).
Я настоятельно рекомендую использовать #P - оно не решает само никаких проблем, но сильно помогает узнать в чем же они.
Также советую в следующий раз прилагать аутпут. (если расширение не поддерживается - переименуйте в txt)

[Phoebe12]

По поводу длины схемы согласна с Вами. Просто с изначальным использованием B3LYP/6-311++G геометрия системы почему-то не смогла оптимизироваться, поэтому я решила хотя бы попробовать пойти путём "от меньшего к большему" с надеждой на оптимизацию.

Большое спасибо, воспользуюсь Вашими советами! Если после этого проблема не решится, обязательно приложу итоговый output.

[Гесс]

кстати не обратил внимания. Тут viewtopic.php?f=71&t=112356&hilit было сломано много копий, но я по прежнему при своем мнении и считаю что в вашей системе, особо при зафиксированных некоторых атомах из рентгена, без использования D3 поправки можно получить ересь. (водородная связь разумеется кое-как описывается и без поправки гримме, но тут слишком много межмолекулярки чтобы совсем ей пренебрегать).

[Phoebe12]

Спасибо, почитаю и попробую вникнуть
У меня тут ещё в ходе изменения input-файла для сканирования один примитивный вопрос назрел... Когда я построила свою исходную систему атомов в GaussView, он выдал, что заряд +1 и это дублет, и я особо не задумалась (а зря) и задала эти характеристики в input файле. Но если в системе 3 молекулы воды и 2 незаряженных полярных аминокислоты, то справедливо ставить 0 и 1 соответственно? А то какая-то ересь получается

[amge]

Обратите внимание на 45-й протончик. Он как раз и есть Н+. Посредине между кислородами, что для такого типично (Zindel cation).

[Phoebe12]

Логично... Тогда у меня путаница в понимании следующего: если брать вообще голую систему (то есть до предварительной оптимизации), когда есть тирозин с группой OH на своём кольце (с которого и будет переноситься протон в ходе сканирования), есть две незаряженные аминокислоты и три молекулы воды, то такая система имеет заряд 0 и мультиплетность 1? Или я ошибаюсь?

[Гесс]

А это вам решать чего вы хотите. Можно незаряженную систему без неспаренных электронов 0 1
А можно забрать у системы один электрон и считать катионрадикал 1 2
Вам скорее всего интересен первый случай. Но за время оптимизации заряд и мультиплетность системы в целом - не менялись.

[Phoebe12]

Последовала Вашему совету и сделала релаксированный скан на B3LYP/6-31+G* с подробным аутпутом (прилагаю ниже инпут в виде .txt, а аутпут по ссылке ниже на Гугл Диск, так как он не захотел заливаться на форум из-за размера файла).

Печаль началась уже как минимум на 3 шаге сканирования, когда геометрия не смогла оптимизироваться, затем то же произошло и на 4-5 шагах, а в конечном итоге появилась Server Error Message 2070. В связи с этим буду рада услышать ответы/советы по поводу нескольких вопросов:

1) В чём может быть проблема с оптимизацией на каждом шаге сканирования в этот раз и как её можно устранить?
2) Какими путями можно и нужно ускорить расчёты (понимаю, что 4 дня - это не месяц и даже не 2 недели, но всё-таки)?
3) Вдогонку ко 2 пункту: есть ли какие-то варианты квантовохимических пакетов, которые считают DFT быстрее Гауссиана?

Ссылка на Гугл Диск, где лежит аутпут: https://drive.google.com/file/d/0B8l7ut ... sp=sharing

[Гесс]

Таак. Ну начнем с конца.
3) Быстрее чем Гауссиан негибридные функционалы ситает Природа (помнится по умолчанию без RI).
Быстрее чем гауссиан гибридные функционалы может считать ОРКА, с RI, просто потому что у гауссиана нет RI для гибридных функционалов.
Ну и негибридные функцционалы можно ускорить с помощью RI в гауссиане.

2) Оно настолько долго потому что оно делает по три сотни шагов оптимизации на точку скана вместо трех десятков. Если вылечить сходимость эти расчеты будут проскакивать за один день. Разве что свопится от нехватки памяти, 1 Гб на 8 ядер это немного.

1) ну совету касательно исходной проблемы вы не последовали. ИМХО там идут какие-то мелкие колебания (на третьем шаге четко просматривается классическая осцилляция в конце). Во первых я по прежнему предлагаю использовать декартовы координаты, это решит проблему если она связана со внутренними. (декартовы координаты не работают в паре со сканом поэтому придется делать отдельные расчеты). Во вторых я по прежнему предлагаю использовать дисперсионную поправку. Ну и в третьих можно попробовать улучшить сетку интегрирования, но это бы я оставил на потом, если первые два совета не сработают.

[Phoebe12]

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

Можете расшифровать, пожалуйста, что имеется ввиду под отдельными расчётами?

[Гесс]

Я имел в виду вместо скана сделать 5 независимых расчетов с констрейнутыми геометриями, но уже не важно. Я это уже проверил и точка 5 не работает весьма сходно тому как не работает у вас. Optimization stopped. -- Number of steps exceeded, NStep= 148
И проблемы со сходимостью на последних шагах. Это надо дальше думать. Не сейчас.

[Phoebe12]

А если попробовать обойтись без Relaxed PES, оптимизировав констрейнутую геометрию отдельно для реагентов и продуктов, а потом запустить QST2? Или там всё может быть ещё хуже?