ValueError in g-tensor calculation

[1234zou]

Dear PySCF developers,
I'm trying to use this repo to perform a UHF g-tensor calculation. But got the following error

Traceback (most recent call last):
  File "/public/home/jxzou/vip_test/Tp2Co_C6H6/Co_quartet/Co_quartet_uhf_nmr.py", line 22, in <module>
    gobj.kernel()
  File "/public/home/jxzou/software/properties/pyscf/prop/gtensor/uhf.py", line 486, in kernel
    gpara = self.para(mo1, self._scf.mo_coeff, self._scf.mo_occ)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/public/home/jxzou/software/properties/pyscf/prop/gtensor/uhf.py", line 513, in para
    self.mo10, self.mo_e10 = self.solve_mo1()
                             ^^^^^^^^^^^^^^^^
  File "/public/home/jxzou/software/properties/pyscf/prop/nmr/uhf.py", line 174, in solve_mo1
    mo10, mo_e10 = ucphf.solve(vind, mo_energy, mo_occ, h1, s1,
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/public/home/jxzou/software/pyscf-2.11.0/pyscf/scf/ucphf.py", line 41, in solve
    return solve_withs1(fvind, mo_energy, mo_occ, h1, s1,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/public/home/jxzou/software/pyscf-2.11.0/pyscf/scf/ucphf.py", line 147, in solve_withs1
    mo1 = lib.krylov(vind_vo, mo1base.reshape(-1, nov),
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/public/home/jxzou/software/pyscf-2.11.0/pyscf/lib/linalg_helper.py", line 1311, in krylov
    axt = aop(x1)
          ^^^^^^^
  File "/public/home/jxzou/software/pyscf-2.11.0/pyscf/scf/ucphf.py", line 137, in vind_vo
    v = fvind(mo1).reshape(nd, nov)
        ^^^^^^^^^^
  File "/public/home/jxzou/software/properties/pyscf/prop/nmr/uhf.py", line 203, in vind
    v1mo = numpy.hstack((numpy.asarray(v1a),
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/public/home/jxzou/software/anaconda3/lib/python3.11/site-packages/numpy/core/shape_base.py", line 359, in hstack
    return _nx.concatenate(arrs, 1, dtype=dtype, casting=casting)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 2, the array at index 0 has size 15 and the array at index 1 has size 12

The input script is

from pyscf import gto, scf, lib
from pyscf.prop import gtensor

lib.num_threads(32)
mol = gto.M()
mol.atom = '''
Co           0.000000000         0.000000000         0.000000000
'''

mol.basis = 'def2-TZVP'
mol.charge = 0
mol.spin = 3
mol.verbose = 4
mol.build(parse_arg=False)

mf = scf.UHF(mol)
mf.max_cycle = 128
mf.max_memory = 64000 # MB
mf.kernel()

gobj = gtensor.uhf.GTensor(mf).set(verbose=4)
gobj.kernel()

Is the problem related to numpy version? Here is what I got from conda list numpy

# packages in environment at /public/home/jxzou/software/anaconda3:
#
# Name                    Version                   Build  Channel
numpy                     1.26.4          py311h08b1b3b_0
numpy-base                1.26.4          py311hf175353_0
numpydoc                  1.5.0           py311h06a4308_0

Python version is 3.11.7.

Thanks a lot,
Jingxiang

[jeanwsr]

This is incorrect since we cannot hstack two array with shape ( .., nao, na/nb) when na not equals nb.

properties/pyscf/prop/nmr/uhf.py

Lines 203 to 204 in 4eee5a4

	v1mo = numpy.hstack((numpy.asarray(v1a),
	numpy.asarray(v1b)))

[jeanwsr]

Could you try if #13 fixes the problem?

[1234zou]

Thanks a lot @jeanwsr , The hstack error is fixed. For the Co atom, the Gaussian result is

 g value of the free electron [g_e]: 0.20023193D+01
 relativistic mass correction [g_RMC]:-0.53351695D-02
 diamagnetic correction to g tensor [g_DC]:
   XX= 0.39904788D-02   YX= 0.30829295D-13   ZX= 0.34218577D-11
   XY= 0.30829295D-13   YY= 0.39904788D-02   ZY= 0.10254579D-11
   XZ= 0.34218577D-11   YZ= 0.10254579D-11   ZZ= 0.39904788D-02
 orbital Zeeman and spin-orbit coupling contribution to g tensor [g_OZ/SOC]:
   XX= 0.10269727D+02   YX=-0.21955117D-08   ZX= 0.79931501D-07
   XY=-0.21955065D-08   YY= 0.10269727D+02   ZY= 0.23499979D-07
   XZ= 0.79931465D-07   YZ= 0.23500006D-07   ZZ= 0.10269727D+02
 g tensor [g = g_e + g_RMC + g_DC + g_OZ/SOC]:
   XX= 0.12270702D+02   YX=-0.21954809D-08   ZX= 0.79934923D-07
   XY=-0.21954757D-08   YY= 0.12270702D+02   ZY= 0.23501004D-07
   XZ= 0.79934886D-07   YZ= 0.23501031D-07   ZZ= 0.12270702D+02
 g shifts relative to the free electron (ppm):
   xx=     10268382.5   yy=     10268382.6   zz=     10268382.7

Note that the off-diagonal elements are not exactly zero since I usually use nosymm in any calculation. If I remove this keyword, the off-diagonal elements become either 0.0 or 1e-33. So I assume that only diagonal elements are useful. The g_OZ/SOC term of Gaussian is vastly different from those of ORCA/PySCF (I think there is probably some bug in this term for Gaussian).

The ORCA result (using ! UHF VeryTightSCF noTRAH ZEFF-SOC) is

The g-matrix:
              2.6825670   -0.0000000    0.0000000
             -0.0000000    2.6825670    0.0000000
              0.0000000    0.0000000    2.6825670

Breakdown of the contributions
 gel          2.0023193    2.0023193    2.0023193
 gRMC        -0.0005935   -0.0005935   -0.0005935
 gDSO(tot)    0.0002645    0.0002645    0.0002645
 gPSO(tot)    0.6805767    0.6805767    0.6805767
             ----------   ----------   ----------
 g(tot)       2.6825670    2.6825670    2.6825670 iso=  2.6825670
 Delta-g      0.6802477    0.6802477    0.6802477 iso=  0.6802477

And the PySCF result from #13 is

dia_soc2e = False
para_soc2e = SSO+SOO
so_eff_charge = False (1e SO effective charge)
RMC = -0.0005927972519884573
GC(1e) [0.00026415 0.00026415 0.00026415]
SOC(1e)/OZ [1.21394779 1.21391927 1.21397072]
SOC(2e)/OZ [-0.52263662 -0.52262429 -0.52264653]
free electron g 2.00231930436182
g-tensor diamagnetic terms [-0.00032865 -0.00032865 -0.00032865]
g-tensor paramagnetic terms [0.69131117 0.69129498 0.69132419]
g-tensor total [2.69330182 2.69328564 2.69331485]
g-shift (ppt) [690.98251522 690.96633231 690.99554154]

It seems that the g-tensor diamagnetic terms of PySCF has an opposite sign. The similar sign problem can also be observed for the doublet Sc atom. But their absolute values are small, I'm not sure whether it is a numerical error problem. Getting consistent g-tensors from three software packages may be beyond the scope of this issue. If you think PySCF g-tensor diamagnetic terms is correct, I think this issue can be closed.