Use of DOMONLY and SAVEDOM in PNO-LCCSD calculations

[Liz D]

As per https://pubs.acs.org/doi/10.1021/acs.jctc.8b01098 (and the molpro documentation on older local methods), the "best practice" for getting a smooth PES over a scan of intermolecular distances is to determine the domains at one fixed separation (e.g., the equilibrium structure) and freeze the domains there for all other distances.

I am trying to do this in Molpro 2019.2, and cannot seem to make the options work.

Specifically, both:

local,domopt=tight
...
{pno-lccsd, implementation=disk
local,domonly=1,savedom=6200.2}

And

local,domopt=tight
...
{pno-lccsd, implementation=disk}

Provide the same output: a full PNO-LCCSD calculation with no hint anywhere in the output that domains are being saved to 6200.2. Molpro does seem to think my local options are valid (it crashes with a rejection if I change instead to some gibberish down there), but... as far as I can tell nothing is happening once those options are parsed. I've tried putting domonly=1 in the first local (nothing happens) and domopt=tight in the pno-lccsd's local (molpro rejects this).

When I flip to df-lccsd, I get the results I expect (domains stored on the record, calculation exits without actually performing a CCSD-level calculation)

What is the proper way to format my file to: 1. calculate only domains and then move on, and 2. save those domains somewhere useful, when using PNO?

[qia...@theochem.uni-stuttgart.de]

At present, the domain save/restart feature in PNO-LCCSD is highly experimental. PNOs are more complicated than the PAOs and we have observed (presumably) numerical issues in certain cases. Nevertheless, these calculations can be performed with

{pno-lccsd(t)-f12, save=5000.2, ...}
...
{pno-lccsd(t)-f12, start=5000.2, loc_method=ov, ...}

In the restart calculation you should see output like "Reading reference orbitals from record  5000.2" indicating that the domain restart is in effect.

[Hans-Joachim Werner]

domonly only works for the old PAO local methods. See Qianli Ma’s reply regarding PNO methods.

---
Prof. Dr. Hans-Joachim Werner
Institut für Theoretische Chemie
Universität Stuttgart
Pfaffenwaldring 55
70569 Stuttgart, Germany
e-mail: wer...@theochem.uni-stuttgart.de

> --
> You received this message because you are subscribed to the Google Groups "molpro-user" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to molpro-user...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/molpro-user/6400f9d1-5ab6-4261-aacf-bd18596d5f3fo%40googlegroups.com.

[Liz D]

Yes, that seems to have worked, thank you!

I now have a different problem - I think I am picking a bad "reference" for me to save domains from. In the case of S66x8, point '0' is close to the equilibrium separation, so it sticks out as an obvious "reference". But when picking a random dimer geometry to start with (e.g., as extracted from MD), where I scan intermolecular distances, I do not know a-priori what the equilibrium distance is. I tried two approaches - picking the MD distance as the reference to calculate domains at, and picking a very large intermolecular distance as the reference to calculate domains at.

Neither produced satisfactory results in my test system, a methanol dimer. Using canonical CCSD(T*)-F12b/aVTZ-F12 as my reference, running PNO-LCCSD(T*)-F12b/aVTZ with no freezing of domains gave me a MAD in post-HF energy of 0.01 kcal/mol across all points in my scan, even close-contact points. The only reason I am hesitant to use these results is that the PNO-LCCSD(T*)-F12b/aVTZ data is noticeably non-smooth, where a change in the assignment of pairs to strong/close/weak (I think?) is causing a discontinuity. When using the MD configuration as the reference, the MAD increases to 0.06 kcal/mol, with errors increasing as one moves away from the MD geometry. When using a reference of a very distant intermolecular distance, the MAD is 6.1 kcal/mol, with the error growing as the molecules come closer together (not surprisingly). 

While I could attempt to merge domains along the scan path, I'm not sure how to do this efficiently given that DOMONLY doesn't work and given that I (think) I'd have to double back through the scan once or twice. Nor am I even sure this is the correct approach to take. I could also try taking the closest-contact as my reference, but that feels like a dangerous game given that "closest" may be very close indeed (and in some systems will involve hydrogen bonding).

Do you have any advice on how to proceed? I can tolerate the bumpiness in the non-frozen approach but I suspect that I can do better, though I do not yet know how.

[qia...@theochem.uni-stuttgart.de]

From our experience the discontinuities in intermolecular potentials are mostly caused by the pair approximations. If your molecules are not very large, a good starting point would be turning off the close pair approximations (thrclose=0) and see if the cost and computed energies are acceptable. The domain freezing approach is implemented mostly for numerical energy derivatives and is not suitable when there is a significant change in the molecular geometry.

Using the pair list (and not the domains) determined from the geometry with shortest intermolecular distance is safe. It is unclear how much this can save, because most pairs would be strong in the "closest" geometry if the monomers are not big. Unfortunately this option is not in the program at present (we can consider adding it in the future).

[Liz D]

Given that I intend to work with larger molecules in the future, setting thrclose=0 feels like a hotfix that wouldn't scale for the future. That said, I *can* tolerate the "bumpiness" in the PES, as long as I can understand where it is coming from and feel comfortable that it's not a sign of me doing something "wrong". Everything you've said is very reasonable, and I don't see any need to add the saving of pairs as you're right; the pair list from the closest would be both computationally expensive (and likely overkill!) for most circumstances.

One feature I'm curious if it's possible to work with: is there any way to perturb the pair list during a calculation? E.g., I'd like to know what *change* in energy I can expect from taking the most distant pair in the 'strong' list and shifting it to the 'close' list. I would like to use this to attempt to bound the precision I expect out of PNO for the systems I'm dealing with (I think) when it comes to these scans. Ideally I'd like to redo as little work as possible when it comes to the quantum calculation (even if it involves a non-self-consistent estimate of the energy with the one pair moved from 'strong' to 'close'). While I *think* I could do this by resetting thrclose and re-running within a run, I don't know what the 'magic' value would be to kick exactly one pair down or up the hierarchy. 

[qia...@theochem.uni-stuttgart.de]

Unfortunately it is difficult to estimate accurately the magnitude of the jump caused by a pair class change. The "first order" change comes from the correlation energy of the pair itself, which we know is around 1e-4 Eh (thrclose), or 0.06 kcal/mol with default settings. But we cannot easily tell how much that changes the energies of other pairs through coupling, which can be more than the energy of the pair itself. On the other hand, the close-pair approximations are much better than neglecting the pairs, and this should compensate for the errors a bit. If the question is on the pairs close to the default cut-off threshold, I think it is fine using an empirical bound of error (e.g. the magnitude of the jumps you see on typical PESs) and that is not expected to change too dramatically in similar systems.

[Xiansheng Gao]

Hello, when I use PNO-CCSD(T)-F12 to calculate an open-shell molecule and attempt to obtain gradient information, I specify {pno-lccsd(t)-f12, save=5000.2} in the input. However, the calculation stops with the following error message. Is this functionality has not yet been implemented?  

 Pair list saved on record          5000.2
 ? Error
 ? restart for open-shell molecules is not implemented
 ? The problem occurs in save_orbdom