F/I-SAPT for (formally) open-shell fragments

Hello All,

After being pointed towards the (relatively) new F/I-SAPT capabilities in PSI4, I am in the process of setting up my first intramolecular interaction energy analysis. I want to check how two parts of a closed-shell molecule interact with each other.

As mentioned in the manual (http://www.psicode.org/psi4manual/master/fisapt.html), the examples there do not actually work in the present version, as they define unphysical closed shell singlets with odd number of electrons. There is an example on ISAPT using three fragments (https://github.com/psi4/psi4/blob/master/samples/isapt1/input.dat) but with perhaps somewhat incomplete documentation.

I would highly appreciate it if someone more experienced could say if the input for two fragments below is properly defined. It is supposed to compute the interaction of the methyl group (the first four atoms) with the rest of the molecule. It does run, but my untrained eye cannot say if the numbers are reasonable (except for the dispersion interaction, which seems to be about twice the DFT-D3 value).

Of course, if it is technically correct, but for some other reason a meaningless or dubious way of computing the interaction in question, let me know :slight_smile:

Grateful for any comments, cheers,
Mikael J.

memory 48 GB

#molecule mol {
#0 1
#C         -2.93552        0.29807        3.11475
#H         -3.04406        1.25120        3.68240
#H         -3.76958        0.26173        2.37890
#H         -3.08368       -0.54537        3.82856
#--
#0 1
#C         -5.47668       -0.32237       -0.07427
#C         -5.66953        0.99702       -0.47472
#C         -4.56353        1.77357       -0.81641
#C         -3.26561        1.26867       -0.69202
#C         -3.05378       -0.03702       -0.22022
#C         -4.18115       -0.82923        0.03478
#H         -6.34162       -0.95836        0.17644
#H         -6.68607        1.42658       -0.52901
#H         -4.72861        2.80792       -1.17141
#H         -4.05009       -1.87755        0.35161
#Cl        -1.95937        2.32476       -1.11505
#C         -1.68141       -0.64013        0.00255
#H         -1.62090       -1.58021       -0.59604
#C         -1.40383       -0.97221        1.48100
#H         -0.86421        0.01704       -0.37504
#C         -1.54739        0.20885        2.45863
#H         -2.06649       -1.80502        1.81356
#H         -0.36243       -1.37522        1.54970
#H         -0.80409        0.09928        3.28655
#H         -1.30357        1.16735        1.94440
#symmetry c1
#no_reorient
#no_com
#}

mol = psi4.core.Molecule.from_arrays(
    elez=[6, 1, 1, 1, 6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 17, 6, 1, 6, 1, 6, 1, 1, 1, 1],
    fragment_separators=[4],
    fix_com=True,
    fix_orientation=True,
    fix_symmetry='c1',
    fragment_multiplicities=[2, 2],
    molecular_charge=0,
    molecular_multiplicity=1,
    geom= [
     -2.93552   ,    0.29807  ,     3.11475,
     -3.04406   ,    1.25120  ,     3.68240,
     -3.76958   ,    0.26173  ,     2.37890,
     -3.08368   ,   -0.54537  ,     3.82856,

     -5.47668   ,   -0.32237  ,    -0.07427,
     -5.66953   ,    0.99702  ,    -0.47472,
     -4.56353   ,    1.77357  ,    -0.81641,
     -3.26561   ,    1.26867  ,    -0.69202,
     -3.05378   ,   -0.03702  ,    -0.22022,
     -4.18115   ,   -0.82923  ,     0.03478,
     -6.34162   ,   -0.95836  ,     0.17644,
     -6.68607   ,    1.42658  ,    -0.52901,
     -4.72861   ,    2.80792  ,    -1.17141,
     -4.05009   ,   -1.87755  ,     0.35161,
     -1.95937   ,    2.32476  ,    -1.11505,
     -1.68141   ,   -0.64013  ,     0.00255,
     -1.62090   ,   -1.58021  ,    -0.59604,
     -1.40383   ,   -0.97221  ,     1.48100,
     -0.86421   ,    0.01704  ,    -0.37504,
     -1.54739   ,    0.20885  ,     2.45863,
     -2.06649   ,   -1.80502  ,     1.81356,
     -0.36243   ,   -1.37522  ,     1.54970,
     -0.80409   ,    0.09928  ,     3.28655,
     -1.30357   ,    1.16735  ,     1.94440])
activate(mol)

set {
basis         jun-cc-pvdz
scf_type df
guess sad
freeze_core true
}

energy('fisapt0')