SAPT(0) Calculation for transition metal interaction

SAPT
1

Dear colleagues, I’m trying the excellent SAPT procedures in Psi4 to evaluate interactions between a ligand and a transition metal complex. As a toy system, I’m using nickelocene:

memory 2 GB

molecule {
-1 1
   C        0.00000        0.00000        2.14426
   C        1.30823        0.00000        1.71870
   C        1.54069        1.13775        1.04373
   H        2.02869       -0.78448        1.89806
   C        0.41917        1.88417        0.99981
   C       -0.60852        1.20803        1.68500
   H       -0.48098       -0.77477        2.72230
   H        0.32449        2.84521        0.51681
   H       -1.62923        1.53050        1.83013
   H        2.48626        1.41305        0.60172
--
1 3
  Ni        0.00000        0.00000        0.00000
   C        0.00000        0.00000       -2.14426
   C       -1.30823        0.00000       -1.71870
   C       -1.54069       -1.13775       -1.04373
   C       -0.41917       -1.88417       -0.99981
   C        0.60852       -1.20803       -1.68500
   H        0.48098        0.77477       -2.72230
   H       -2.02869        0.78448       -1.89806
   H       -2.48626       -1.41305       -0.60172
   H       -0.32449       -2.84521       -0.51681
   H        1.62923       -1.53050       -1.83013

symmetry c1
no_reorient
}
set {
	reference uhf
	basis cc-pVDZ
	maxiter 2000
}
energy('sapt0')

The SCF part is converging nicely, but at the CPKS Iterations the calculation dies:

==> CPKS Iterations <==

    Maxiter     =        2000
    Convergence =   1.000E-09

    -----------------------------------------
    Iter   Monomer A    Monomer B    Time [s]
    -----------------------------------------
    1      1.454E-01    4.016E-01           6
    2      5.764E-02    1.471E-01          12
    3      9.373E-03    5.043E-02          18
    4      4.788E-03    2.959E-02          23
    5      1.404E-03    1.228E-02          29
    6      2.083E-04    6.786E-03          35
    7      1.131E-04    3.878E-03          41
    8      2.803E-05    3.730E-03          46
    9      1.097E-05    4.778E-03          52
    10     2.086E-06    4.808E-03          58
    11     1.188E-06    5.031E-03          63
    12     3.698E-07    5.783E-03          68
    13     7.100E-08    5.159E-03          73
    14     5.979E-08    1.679E-02          78

Traceback (most recent call last):
  File "/home/henrique/bin/anaconda3/envs/psi4/bin/psi4", line 287, in <module>
    exec(content)
  File "<string>", line 52, in <module>
  File "/home/henrique/bin/anaconda3/envs/psi4/lib//python3.7/site-packages/psi4/driver/driver.py", line 556, in energy
    wfn = procedures['energy'][lowername](lowername, molecule=molecule, **kwargs)
  File "/home/henrique/bin/anaconda3/envs/psi4/lib//python3.7/site-packages/psi4/driver/procrouting/proc.py", line 3404, in run_sapt
    e_sapt = core.sapt(dimer_wfn, monomerA_wfn, monomerB_wfn)

RuntimeError: 
Fatal Error: Monomer B: A Matrix is not SPD
Error occurred in file: /scratch/psilocaluser/conda-builds/psi4-multiout_1557940846948/work/psi4/src/psi4/libsapt_solver/usapt0.cc on line: 1465
The most recent 5 function calls were:




Printing out the relevant lines from the Psithon --> Python processed input file:
    """)
    core.IO.set_default_namespace("")
    core.set_global_option("REFERENCE", "uhf")
    core.set_global_option("BASIS", "cc-pVDZ")
    core.set_global_option("MAXITER", 2000)
--> energy('sapt0')

    Psi4 stopped on: Wednesday, 05 August 2020 04:40PM
    Psi4 wall time for execution: 0:06:19.30

*** Psi4 encountered an error. Buy a developer more coffee!
*** Resources and help at github.com/psi4/psi4.

I’m using Psi4 1.3.2 release from Anaconda.

2

That error message means that you haven’t landed on the correct ground state for the second monomer. (More precisely, you’re not at a local minimum in orbital-space.)

The normal tricks I’d use to improve the SCF convergence aren’t working. If I turn on SOSCF, the convergence is preposteorusly slow. If I don’t turn on SOSCF but do use stability analysis, the iterations stall out before 1e-7 RMS density error for some reason.

The best I can recommend is that you look http://psicode.org/psi4manual/1.3.2/scf.html, including the “Keywords” section, and start trying things. Others may have more helpful suggestions.

3

You may want to look at alternative choices for partitioning the molecule into fragments. This is stretching my knowledge of inorganic chemistry, but I think the Ni electrons live in d orbitals that interact with both of the rings. Thus, it’s not appropriate to associate the Ni with one ring but not the other, as you’ve done in your input file. I think this is why you are having trouble converging the SCF procedure for the monomer.

I think your best option is intramolecular SAPT, which will use a more rigorous approach to localize orbitals to particular functional groups. You should use the input format described here in which you specify the charge & multiplicity of the entire molecule first and then specify charge & multiplicity of each fragment.

You will still have to make a choice about how to partition the molecule into fragments. It is usually recommended to avoid cutting fragments across bonds with pi character and to avoid having a single atom by itself, but I don’t know if transition metals are an exception. Maybe @dasirianni can weigh in.

In any case, after running fsapt.py, you should check the “Orbital Check” section of the output. If any value is bigger than 0.1, then you’ve made a poor choice in partitioning the molecule.

4

Hello, and sorry for my huge delay. Thank you for your insightful replies, it really helped me to achieve some progress but at the end it looks like I’m not going to be able to use Psi4 for my work as the sapt0 methodology failed to converge for all of my open-shell cases.
fisapt0 says that it is defined only for closed shell cases, so I wasn’t able to use it either.

Who knows in a future update? For now all I can do is envy some of the Psi4 features and hope it works for my systems soon. :smiley: