ADMP в Gaussian 03

[fullerene]

Коллеги!

Объясните неспециалисту, почему при расчете маленькой системы (4-20 атомов) методом ADMP в Gaussian03 время расчета существенно увеличивается, если задать трансляционные векторы (TV). Размер бокса, задаваемого векторами, сопоставим или чуть больше размера системы, рассчитанного по координатам крайних атомов. Вот эти параметры не задаются: http://www.gaussian.com/g_tech/g_ur/k_pbc.htm

[Droog_Andrey]

http://python.rice.edu/~guscus/gaussian ... guide.html

Попробуйте уменьшить базис (например, убрав диффузные функции), уменьшить диапазон поиска ghost-атомов (например, задав PBC(CellRange=50)).

[fullerene]

Вот такая конкретная задача:
# admp=(stepsize=2500,maxpoints=20) ub3lyp/3-21g temperature=2000 guess=(read,fock)
Система: NO + O₂, энергия минимизирована.

Имеем:
PBC vector 1 X= 7.5589 Y= 0.0000 Z= 0.0000
PBC vector 2 X= 0.0000 Y= 7.5589 Z= 0.0000
PBC vector 3 X= 0.0000 Y= 0.0000 Z= 7.5589
Recp vector 1 X= 0.1323 Y= 0.0000 Z= 0.0000
Recp vector 2 X= 0.0000 Y= 0.1323 Z= 0.0000
Recp vector 3 X= 0.0000 Y= 0.0000 Z= 0.1323
One-electron integrals computed using PRISM.
NBasis= 36 RedAO= T NBF= 36
NBsUse= 36 1.00D-06 NBFU= 36
RepCel: MaxNCR= 257 NClRep= 257 NMtPBC= 1.
Initial guess read from the checkpoint file:
C:\G03W\MySystems\NO2_3.chk
<S**2> of initial guess= 1.6303
ReFckM: NMPBCI= 1 NMPBCO= 1 Same=T
Requested convergence on RMS density matrix=1.00D-07 within 128 cycles.
Requested convergence on MAX density matrix=1.00D-05.
Requested convergence on energy=1.00D-05.
No special actions if energy rises.
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
Rare condition: small coef for last iteration: 0.000D+00
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
EnCoef did 3 forward-backward iterations
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
EnCoef did 100 forward-backward iterations
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
EnCoef did 100 forward-backward iterations
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
EnCoef did 2 forward-backward iterations
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F
FoFJK: IHMeth= 1 ICntrl= 0 DoSepK=T KAlg= 1 I1Cent= 0 FoldK=F

Не пониманию смысл выделенного жирным.

[Droog_Andrey]

SCF не сходится. Возможно, проблема в мультиплетности.

А зачем здесь вообще PBC?

[fullerene]

SCF не сходится

А без TV сходится...

А зачем здесь вообще PBC?

Здесь - низачем. Просто маленькая система для разбора проблемы с большей рабочей системой.

[Droog_Andrey]

А Вы энергию без TV минимизировали?

Впрочем, если Вам так уж понадобились PBC, то нужно для начала учесть пожелания системы по поводу параметров ячейки, Z и т.п.

Маловероятно, что NO...O₂ образует ферромагнитный кристалл, да ещё при Т = 2000К

А в чём состоит проблема с большой системой?

[Darth Vasya]

Коллеги!
время расчета существенно увеличивается, если задать трансляционные векторы (TV) ... Вот эти параметры не задаются: http://www.gaussian.com/g_tech/g_ur/k_pbc.htm

Дык теорема Блоха ж! Гауссиан по умолчанию берёт количество k-точек с потолка - видимо, чтоб заведомо не недобдеть. Задайте NKpoint, соответствующее вашей системе: если диэлектрик, то порядка десятка, если металл - то побольше. Ну и увеличивайте постепенно до достижения сходимости.

Ну и непрошенный, но намного более уместный совет: оставьте вы эту идею считать периодические системы гауссианом, это ни к чему путному не приведёт. Возьмите подходящий софт. Если обязательно хотите пропагацию волновых функций, то смотрите в сторону CPMD или Quantum ESPRESSO (или CP2k). Но это только если вы уверены, что приближение Борна-Оппенгеймера вам почему-то не подходит. В противном случае - вообще любой пакет для периодических систем подойдёт. Если обязательно нужны гауссовы функции, то Crystal (или CP2k, опять же), если подойдут и плоские волны - могу посоветовать ABINIT, он умеет работать с projector-augmented waves.

Вообще, стремление использовать гауссиан для всего на свете - типичное проявление "закона инструмента" (if all you have is a hammer, everything looks like a nail)

[Darth Vasya]

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

[fullerene]

Ну и непрошенный, но намного более уместный совет: оставьте вы эту идею считать периодические системы гауссианом, это ни к чему путному не приведёт.

Может быть...

Возьмите подходящий софт. Если обязательно хотите пропагацию волновых функций, то смотрите в сторону CPMD или Quantum ESPRESSO (или CP2k). Но это только если вы уверены, что приближение Борна-Оппенгеймера вам почему-то не подходит.

Пользуясь случаем... В CPMD как-то можно выводить заряды на атомах каждые N-шагов в процессе расчета траектории? Или только скриптом с поочередным вызовом & CPMD PROPERTIES || MOLECULAR DYNAMICS CP &END..?

[Darth Vasya]

Что касается CPMD, то с такими подробностями я вряд ли смогу помочь: мне с ним плотно работать не доводилось - так, потыкал когда-то и забыл. Впрочем, для таких вопросов есть официальный CPMD-спамлист Более того, есть все шансы, что этот вопрос там уже обсуждался: http://www.cpmd.org/cpmd_mailinglist.html

По идее, должно работать прямо так:

&CPMD
PROPERTIES
MOLECULAR DYNAMICS
...
&END

...
&PROP
CHARGES
&END

Так должно выводить каждые N=1 шагов, дальше тривиально

И ещё раз разрешите полюбопытствовать: чем именно вас товарищи Борн с Оппенгеймером не удовлетворяют?

[fullerene]

PROPERTIES
MOLECULAR DYNAMICS

В CPMD эти директивы взаимоисключающие. И это странно.

И ещё раз разрешите полюбопытствовать: чем именно вас товарищи Борн с Оппенгеймером не удовлетворяют?

Да вообщем-то удовлетворяют не меньше, чем Кар с Парринелло. Просто надо же все "потыкать". В Гауссиане с варьированием NKPOINT не получилось - пока бросил. Рабочая система простая - ионная(ые) пара(ы) [LI+][BF4-] в метаноле.

Это пропаганда ABINIT?

[Darth Vasya]

Нет, не пропаганда: в CPMD тоже есть Борн-Оппенгеймер, да и в Гауссиане, если на то пошло Есть ещё ESPRESSO и CP2k, я бы просто интереса ради последнее попробовал

Что касается формализма CP, конкретно с точки зрения плоских волн меня немножко напрягает, что количество базисных функций фиксированное, т.е. нельзя варьировать параметры ячейки (а точнее, можно, но возникает систематическая ошибка).

Вот ещё любопытный метод, кстати: http://link.aip.org/link/?JCPSA6/130/224106/1

[fullerene]

Нет, не пропаганда: в CPMD тоже есть Борн-Оппенгеймер, да и в Гауссиане, если на то пошло Есть ещё ESPRESSO и CP2k, я бы просто интереса ради последнее попробовал

Я постараюсь...

Вот ещё любопытный метод, кстати: http://link.aip.org/link/?JCPSA6/130/224106/1

Видел, но сказать совершенно ничего не могу по этому поводу.

[Darth Vasya]

PROPERTIES
MOLECULAR DYNAMICS

В CPMD эти директивы взаимоисключающие. И это странно.

Это чрезвычайно странно: я совершнно уверен, что в CPMD можно было "на лету" считать центры функций Ванье - это тоже PROP, вроде бы...

[fullerene]

PROPERTIES
MOLECULAR DYNAMICS

В CPMD эти директивы взаимоисключающие. И это странно.

Это чрезвычайно странно: я совершнно уверен, что в CPMD можно было "на лету" считать центры функций Ванье - это тоже PROP, вроде бы...

Копаюсь. Просто DIPOLE MOMENT точно не работает налету. Про ф-ии Wannier еще посмотрю в мануале.

[alxyppv]

Копаюсь. Просто DIPOLE MOMENT точно не работает налету. Про ф-ии Wannier еще посмотрю в мануале.

эээ... а как на счет опции "DIPOLE DYNAMICS"?

[alxyppv]

MOLECULAR DYNAMICS [CP, BO, PT, CLASSICAL, FILE [XYZ, NSKIP=N, NSAMPLE=M]]
Section: &CPMD

Perform a molecular dynamics (MD) run. CP stands for a Car-Parrinello type MD.
With the option BO a Born-Oppenheimer MD is performed where the wavefunction
is reconverged after each MD-step. CLASSICAL means that a MD that includes
classical atoms is performed.

If FILE is set, then the trajectory is reread from a le instead of being calculated. This
is useful for performing analysis on a previous trajectory. Can be used in conjunction
with the standard MD options like DIPOLE DYNAMICS and WANNIER; some other
features like LINEAR RESPONSE are also enabled. The trajectory is read from
a le named TRAJSAVED (usually a copy of a previous TRAJECTORY le), or
TRAJSAVED.xyz if XYZ is set. NSKIP and NSAMPLE control the selection of
frames read: the frame read at step ISTEP is NSKIP+ISTEP*NSAMPLE.
Default is CP.

[Darth Vasya]

В конце концов, рассчитайте по готовой траектории, в самом деле Прирост затрат не более, чем в два раза, а если не каждый шаг делать - просто ничтожный по сравнению с самой динамикой...

[fullerene]

MOLECULAR DYNAMICS [CP, BO, PT, CLASSICAL, FILE [XYZ, NSKIP=N, NSAMPLE=M]]
Section: &CPMD

Perform a molecular dynamics (MD) run. CP stands for a Car-Parrinello type MD.
With the option BO a Born-Oppenheimer MD is performed where the wavefunction
is reconverged after each MD-step. CLASSICAL means that a MD that includes
classical atoms is performed.

If FILE is set, then the trajectory is reread from a le instead of being calculated. This
is useful for performing analysis on a previous trajectory. Can be used in conjunction
with the standard MD options like DIPOLE DYNAMICS and WANNIER; some other
features like LINEAR RESPONSE are also enabled. The trajectory is read from
a le named TRAJSAVED (usually a copy of a previous TRAJECTORY le), or
TRAJSAVED.xyz if XYZ is set. NSKIP and NSAMPLE control the selection of
frames read: the frame read at step ISTEP is NSKIP+ISTEP*NSAMPLE.
Default is CP.

Спасибо. Разобрался.

[fullerene]

Разобраться-то с мануалом я разобрался, а на практике так и не могу найти, где CPMD выводит заряды на каждом атоме (если он их вобще выводит).

Со строкой DIPOLE DYNAMICS WANNIER создаются куча файлов WC_SPREAD, WC_QUAD, WANNIER_CENTER, WANNIER_HAM, WANNIER_DOS, DIPOLE, но где искать собственно заряды я так и не соображу.