H2+ PDFT Inversion

import numpy as np
import matplotlib.pyplot as plt

from CADMium import Pssolver
from CADMium import Psgrid
from CADMium import Partition
from CADMium import Inverter

a = 2.0/2
#Nuclear charge for fragments A and B
Za, Zb = 1,1
#Set polarization 1-Unpolarized, 2-Polarized|
pol = 2
#Fragment a electrons [alpha, beta]
Nmo_a = [[1  ,0]] #Number of molecular orbitals to calculate
N_a   = [[0.5,0]]
#Ensemble mix
nu_a = 1
#Fragment b electrons
Nmo_b = [[1  ,0]]
N_b   = [[0.5,0]]
#Ensemble mix
nu_b = 1

#Molecular elctron configuration
Nmo_m = [[1,1]]
N_m   = [[1,1]]

#Set up grid
NP = 7
NM = [4,4]
L = np.arccosh(10/a)
loc = np.array(range(-4,5)) #Stencil outline

grid = Psgrid(NP, NM, a, L, loc)
grid.initialize()

part = Partition(grid, Za, Zb, pol, Nmo_a, N_a, nu_a, Nmo_b, N_b, nu_b, {  "ab_sym"            : True,
                                                                           "ens_spin_sym"      : False,
                                                                           "kinetic_part_type" : "libxcke",
                                                                           "k_family"          : "gga",
                                                                           "ke_func_id"        : 500,
                                                                           "interaction_type"  : "ni",
                                                                           "fractional"        : True,
                                                                            })

#Setup inverter object
mol_solver = Pssolver(grid, Nmo_m, N_m)
part.inverter = Inverter(grid, mol_solver, {  "AB_SYM"         : True,
                                              "ENS_SPIN_SYM"   : False,
                                              "use_iterative"  : False,
                                              "invert_type"    : "wuyang",
                                              "Tol_lin_solver" : 1e-3,
                                              "disp"           : True,
                                            })

Warning: If len(KS) > 1 Has not been migrated from matlab

### Perform Isolated Fragment Calculation.
part.optPartition.isolated = True
part.scf({"disp"     : True,
          "e_tol"    : 1e-7})

                  Total Energy

                __________________

Iteration         A            B              res

_______________________________________________________

    1           -0.25000     -0.25000       1.000e+00
    2           -0.25000     -0.25000       8.199e-10

### Perform PDFT Calculation.

part.optPartition.isolated = False
part.scf({'iterative' : False,
          'disp'      : True,
           'continuing' : True})

                  Total Energy

                __________________

Iteration         A            B              res

_______________________________________________________

    1           -0.20949     -0.20949       1.000e+00
    2           -0.21997     -0.21997       7.238e-02
    3           -0.22295     -0.22295       1.886e-02
    4           -0.22339     -0.22339       2.846e-03
    5           -0.22347     -0.22347       5.259e-04
    6           -0.22348     -0.22348       9.523e-05
    7           -0.22349     -0.22349       1.723e-05
    8           -0.22349     -0.22349       3.114e-06
    9           -0.22349     -0.22349       5.623e-07
   10           -0.22349     -0.22349       1.015e-07
   11           -0.22349     -0.22349       1.831e-08