Treating Core Electrons as Valence in PNO-LCCSD(T*)-F12b

25 views
Skip to first unread message

liz.deco...@gmail.com

unread,
Dec 1, 2021, 3:00:03 PM12/1/21
to molpro-user
Molpro version 2021.2

I've been relying on PNO as my main CC-like workhorse for a while now, and generally have had very solid results. Which was why I was surprised that my benchmarking on water/Na+ came out so bad (dimer interaction energies off by >0.4 kcal/mol, when usually I see errors like 0.02 kcal/mol). I'll note that water/Cl- works just fine, and water/NH4+ works just fine -- it appears to be something regarding the need to deal with core electrons that's causing PNO to suddenly behave very, very poorly.

What do I mean by PNO? Specifically I mean PNO-LCCSD(T*)-F12b. The explicit correlation is important, as that seems to be where things fail.

For the same geometry -- Na+ and O, in the avqz-f12,na=cvqz-f12 basis, I see that the PNO-LCCSD(T) - CCSD(T) energy difference for the dimer system is -0.00008081 AU. That's a great match, well within what I expect from extensive use of the theory.

But if I add in F12b correlation, the energy difference jumps to -0.00279642 AU. That's too large for us to work with, and an order of magnitude worse than I've seen for any other system I've simulated (including the somewhat comparable water/NH4+). The difference is entirely in the CCSD: PNO-LCCSD-F12b - CCSD-F12b is -0.00283270 AU, while the difference in triples is a mere 0.00003628 AU.

I've set MODOM=0 explicitly in the mpr just in case that wasn't defaulting right, but it doesn't fix the problem. 

If I run the system with Na+ as a dummy, my energy differences become sensible:
 PNO-LCCSD - CCSD = -0.00012698 AU
 PNO-LCCSD-F12b - CCSD-F12b = -0.00044147 AU

If I run the system with water as a dummy, the problem returns:
PNO-LCCSD - CCSD = 0.00014674 AU
PNO-LCCSD-F12b - CCSD-F12b = -0.00170114

What other settings should I be fiddling with? I think I'm doing something wrong with respect to the need to treat the core electrons in Na+ as valence, but can't figure out how (MODOMC seems to change nothing; I've set it to 0 in both my CCSD-F12b and PNO-LCCSD-F12b calculations).

Input file below:
*** MOLPRO INPUT
memory,692,M,GA=0,M
local,domopt=tight
basis=avqz-f12,na=cvqz-f12
symmetry,nosym
orient,noorient
geomtyp=xyz
geometry={
1
Cartesian coordinates in Angstrom units
Na4 0.00000003 0.00000000 -1.11508555
}
core,1
basis=avtz-f12,na=cvtz-f12
set,charge=1
{hf,accu=15,maxit=50
start,atdens
save,2099.2}
basis=avqz-f12,na=cvqz-f12
set,charge=1
{hf,accu=15,maxit=50
start,2099.2
save,2100.2}
{ccsd(t)}
cc_cc=ENERGC
cc_trip=ENERGT(1) - ENERGC
cc_e=ENERGY
{pno-lccsd(t)}
pno_cc=ENERGC
pno_trip=ENERGT
pno_e=ENERGY
{ccsd(t)-f12b,MODOMC=0}
cc_f12_cc=ENERGC
cc_f12_trip=ENERGT(1) - ENERGC
cc_f12_e=ENERGY
{pno-lccsd(t)-f12b,MODOMC=0}
pno_f12_cc=ENERGC
pno_f12_trip=ENERGT
pno_f12_e=ENERGY

SHOW,pno_cc-cc_cc
SHOW,pno_trip-cc_trip
SHOW,pno_e-cc_e

SHOW,pno_f12_cc-cc_f12_cc
SHOW,pno_f12_trip-cc_f12_trip
SHOW,pno_f12_e-cc_f12_e


Relevant results:
 SETTING PNO_F12_CC     =      -161.98206171  AU                              
 SETTING PNO_F12_TRIP   =        -0.00299633  AU                              
 SETTING PNO_F12_E      =      -161.98505805  AU                              
 PNO_CC-CC_CC     =         0.00014276 AU
 PNO_TRIP-CC_TRIP =         0.00000169 AU
 PNO_E-CC_E       =         0.00014445 AU
 PNO_F12_CC-CC_F12_CC     =        -0.00207440 AU
 PNO_F12_TRIP-CC_F12_TRIP =         0.00001588 AU
 PNO_F12_E-CC_F12_E       =        -0.00205852 AU

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

unread,
Dec 7, 2021, 12:41:20 PM12/7/21
to molpro-user

 Correlated core orbitals are surprisingly difficult for PNO methods. Behind the scenes Molpro does the following to reduce the error:

  - Use tighter PNO thresholds for pairs involving correlated core orbitals.

  - Core orbitals are localized separately from the valence orbitals.

To make that work it is necessary to specify the core orbitals in the command for correlation calculations, e.g.,

  {df-cabs; core, mixed}

  {pno-lccsd(t)-f12; core, mixed}

"core, miexed" is equivalent to "core,1" for Na. See

  https://www.molpro.net/develop/manual/doku.php?id=general_program_structure#defining_orbital_subspaces_occ_closed_core_frozen

for details. In the PNO calculation something like "4 outer core orbitals are localized separately" should be shown in the output file.

"core,1" as a global option as in your input has a slightly different meaning. With that option the PNO program will not consider the 2s2p orbitals of Na as regular valence orbitals (and hence the larger errors). In the next major release we will make PNO-LCCSD recognize "core,mixed" (but not "core, 1") given as a global option.

There was a bug regarding to the recognition of core orbitals in calculations with the dummy atoms, which have just been fixed in the release branch. With that the local errors of Na+...H2O with aug-cc-pwCVTZ basis should amount to around 0.2 kcal/mol (domopt=tight). There are several approaches that could potentially reduce this error significantly. We are currently investigating on them.

The default value of `modomc` is 0 since Molpro 2020.2 so that `modomc=0` has no effect. Also note that, by default, canonical CCSD(T)-F12 does not use the same F12 Ansatz as used in PNO-LCCSD(T)-F12. For the comparison of correlation energies the canonical calculations are best done with

  {CCSD(T)-F12, ansatz=3*A(FIX,NOX)}

liz.deco...@gmail.com

unread,
Dec 7, 2021, 1:35:01 PM12/7/21
to molpro-user
I think I am not understanding something. I have changed my mpr file such that:
`core,1` --> `core,mixed` in my global input

AND

added `{df-cabs;core,mixed}` after my last HF call

AND

added `;core,mixed` to the end of all my PNO (and CC!) commands.

To the MPR file listed above. I get the same energies back out as without these edits, at least to within the precision the printing gives me. This is for a lone Na+, which should have no dummy atoms...

Should I expect this to Work Right with 2021.2? 

Reply all
Reply to author
Forward
0 new messages