Passing reference wavefunctions to SAPT

1

Dear PSI4 community,

I am trying to figure out:
How do I successfully refer previously computed monomer wavefunctions to a SAPT computation?
How can I specify the level of SAPT with psi4.sapt() ?

Basically it should be something like in this example from the manual:
http://psicode.org/psi4manual/master/sapt.html#advanced-example
Unfortunately, whenever I run the code in that example (only copy/pasted from there, no changes), the monomer SCF computations run ok but an error occurs at the last command, the SAPT stage.

Could not find requested basisset (DF_BASIS_SAPT).
Traceback (most recent call last):
  File "/home/tobias/psi4conda/bin/psi4", line 248, in <module>
    exec(content)
  File "<string>", line 69, in <module>

RuntimeError: 
Fatal Error: Wavefunction::get_basisset: Requested basis set (DF_BASIS_SAPT) was not set!

Error occurred in file: /scratch/psilocaluser/conda-builds/psi4_1495011512596/work/psi4/src/psi4/libmints/wavefunction.cc on line: 422
The most recent 5 function calls were:

psi::PsiException::PsiException(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, char const*, int)
psi::Wavefunction::get_basisset(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >)
psi::sapt::USAPT0::USAPT0(std::shared_ptr<psi::Wavefunction>, std::shared_ptr<psi::Wavefunction>, std::shared_ptr<psi::Wavefunction>, psi::Options&, std::shared_ptr<psi::PSIO>)
psi::sapt::sapt(std::shared_ptr<psi::Wavefunction>, std::shared_ptr<psi::Wavefunction>, std::shared_ptr<psi::Wavefunction>, psi::Options&)
py_psi_sapt(std::shared_ptr<psi::Wavefunction>, std::shared_ptr<psi::Wavefunction>, std::shared_ptr<psi::Wavefunction>)

The same happens when I reduce the commands after molecule specification and variables set to (which is basically at the heart of my first question above):

Edim, wfn_dimer = energy('scf',molecule=dimer,return_wfn=True)

monomerA = dimer.extract_subsets(1,2)
EmonA, wfn_monA = energy('scf',molecule=monomerA,return_wfn=True)

monomerB = dimer.extract_subsets(2,1)
EmonB, wfn_monB = energy('scf',molecule=monomerB,return_wfn=True)

psi4.sapt(wfn_dimer,wfn_monA,wfn_monB)

While the error states the basisset is not set, this is in fact the case (I thought? set {df_basis_sapt cc-pVDZ-ri} clearly is present.).I have successfully run more standard SAPT computations on my system, so I guess this is not a local issue.

A problem with this example was pointed out earlier, but no solution was offered (Open-shell SAPT0 computation).

I would greatly appreciate any help in this regard.

Best

2

You are directly calling the C++ layer which is not supported in general. As mentioned you do need to set the basis on the dimer wavefunction object so that it can be used during the SAPT computation. You can find the full SAPT driver here and you can copy and modify as needed.

3

Please let me repeat that this is an example from the manual and it does not run. I cited it here because this is the only reference I could find showing that passing reference wavefunctions to SAPT should be possible, which would be very useful for me. If there is a way to achieve this other than “directly calling the C++ layer” I will gladly use it. I am not a programmer so unfortunately I cannot readily “copy and modify as needed” the SAPT driver with no further help, which would be appreciated.

4

Ah, I see. Please use the following for now. Thank you for bring up this issue, we will post a fix to the documentation:

molecule {
     0 2
     O 0.000000  0.000000  0.000000
     O 0.000000  2.503900  0.000000
     H 0.000000 -0.424700 -1.839500
     --
     0 1
     O 0.000000  0.000000  6.000000
     H 0.000000  1.431500  4.890600
     H 0.000000 -1.431500  4.890600
     units bohr
     symmetry c1
     no_reorient
     no_com
}

dimer = psi4.get_active_molecule()

set {
reference uhf
scf_type     df
basis         cc-pVDZ
df_basis_sapt cc-pVDZ-ri
guess sad
}

dimer = psi4.get_active_molecule()

set df_ints_io save
psi4.IO.set_default_namespace('dimer')
Edim, wfn_dimer = energy('scf',molecule=dimer,return_wfn=True)
set df_ints_io load

monomerA = dimer.extract_subsets(1,2)
psi4.IO.change_file_namespace(97, 'dimer', 'monomerA')
psi4.IO.set_default_namespace('monomerA')
set {
stability_analysis follow
}
EmonA, wfn_monA = energy('scf',molecule=monomerA,return_wfn=True)

monomerB = dimer.extract_subsets(2,1)
psi4.IO.change_file_namespace(97, 'monomerA', 'monomerB')
psi4.IO.set_default_namespace('monomerB')
set {
stability_analysis none
}
EmonB, wfn_monB = energy('scf',molecule=monomerB,return_wfn=True)

psi4.IO.change_file_namespace(97, 'monomerB', 'dimer')
psi4.IO.set_default_namespace('dimer')

aux_basis = psi4.core.BasisSet.build(wfn_dimer.molecule(), "DF_BASIS_SAPT",
                                psi4.core.get_global_option("DF_BASIS_SAPT"),
                                "RIFIT", psi4.core.get_global_option("BASIS"))
wfn_dimer.set_basisset("DF_BASIS_SAPT", aux_basis)
wfn_dimer.set_basisset("DF_BASIS_ELST", aux_basis)

psi4.sapt(wfn_dimer,wfn_monA,wfn_monB)
5

Thank you for your help!