I am wondering if the following script could work.
(1) change basis set add u like:
NickelBasis = BasisSet(
element=PeriodicTable.Nickel,
orbitals=[nickel_3s, nickel_3p, nickel_4s, nickel_3d, nickel_3p_0, nickel_3d_0, nickel_4p],
occupations=[2.0, 6.0, 0.0, 8.0, 1.0, 1.0, 0.0],
hubbard_u=[0.0, 0.0, 0.0, 4.6, 0.0, 4.6, 0.0]*eV,
filling_method=Anisotropic,
onsite_spin_orbit_split=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]*eV,
pseudopotential=NormConservingPseudoPotential("normconserving/sg15/gga/28_Ni.so.upf"),
)
(2) exchange_correlation = ExchangeCorrelation(
exchange=PerdewBurkeErnzerhofExchange,
correlation=PerdewBurkeErnzerhofCorrelation,
hubbard_term=DualShell,
number_of_spins=4,
spin_orbit=True,
)