CCSD(T) 1 particle density matrix

charlotte · May 6, 2026, 8:41pm

Hi,

I try to extract the CCSD(T) 1 particle density matrix for He from psi4.
The following is the code I use, and I get the output showing that the program tried to divide a number by a zero integer. Can anyone help me with this?

import psi4
import numpy as np

# Define the molecule
mol = psi4.geometry("""
0 1
He 0.0 0.0 0.0
symmetry c1
""")

# Set calculation options
psi4.set_options({
    "basis": "STO-3G"
})

# Run the CCSD(T) gradient to trigger the lambda equations and return the wavefunction
_, wfn = psi4.gradient('CCSD(T)', molecule=mol, return_wfn=True)

Da = np.asarray(wfn.Da())
# For RHF-like references, Db may be identical or available; handle both cases:
try:
    Db = np.asarray(wfn.Db())
    P = Da + Db
except Exception:
    P = 2.0 * Da  # common RHF convention if only Da is populated

         ----------------------------------------------------------
                                   FINDIF
                     R. A. King and Jonathon Misiewicz
         ----------------------------------------------------------

  Using finite-differences of energies to determine gradients.
    Generating geometries for use with 3-point formula.
    Displacement size will be 5.00e-03.
    Number of atoms is 1.
    Number of symmetric SALCs is 0.
    Translations projected? 1. Rotations projected? 1.
    Number of geometries (including reference) is 1.
 1 displacements needed ...

  //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>//
  //    FiniteDifference Computations  //
  //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<//

[atl1-1-02-020-33-2:3934722:0:3934722] Caught signal 8 (Floating point exception: integer divide by zero)
==== backtrace (tid:3934722) ====
 0 0x000000000006141c ucs_callbackq_cleanup()  ???:0
 1 0x00000000000616da ucs_callbackq_cleanup()  ???:0
 2 0x000000000003ebf0 __GI___sigaction()  :0
 3 0x0000000000b43927 std::_Hashtable<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::__detail::_Identity, std::equal_to<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::hash<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, std::__detail::_Prime_rehash_policy, std::__detail::_Hashtable_traits<true, true, true> >::_Hashtable<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const*>()  ???:0
 4 0x0000000000b2a681 std::_Hashtable<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::__detail::_Identity, std::equal_to<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::hash<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::__detail::_Mod_range_hashing, std::__detail::_Default_ranged_hash, std::__detail::_Prime_rehash_policy, std::__detail::_Hashtable_traits<true, true, true> >::_Hashtable<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const*>()  ???:0
 5 0x00000000004f9f4a psi::get_block<psi::IntVector>()  ???:0
 6 0x000000000050f3a8 PyInit_core()  ???:0
 7 0x00000000004193a7 psi::get_block<psi::IntVector>()  ???:0
 8 0x00000000001fcdd6 cfunction_call()  /usr/local/src/conda/python-3.11.14/Objects/methodobject.c:542
 9 0x00000000001fcdd6 _Py_CheckFunctionResult()  /usr/local/src/conda/python-3.11.14/Objects/call.c:39
10 0x00000000001fcdd6 cfunction_call()  /usr/local/src/conda/python-3.11.14/Objects/methodobject.c:555
11 0x00000000001d903b _PyObject_MakeTpCall()  /usr/local/src/conda/python-3.11.14/Objects/call.c:214
12 0x00000000001d903b _PyObject_MakeTpCall()  /usr/local/src/conda/python-3.11.14/Objects/call.c:216
13 0x00000000001e73e5 _PyEval_EvalFrameDefault()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:4769
14 0x000000000020d055 _PyEval_EvalFrame()  /usr/local/src/conda/python-3.11.14/Include/internal/pycore_ceval.h:73
15 0x000000000020d055 _PyEvalFrameClearAndPop()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:6406
16 0x000000000020d055 _PyEval_Vector()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:6439
17 0x000000000020d055 _PyFunction_Vectorcall()  /usr/local/src/conda/python-3.11.14/Objects/call.c:393
18 0x00000000002173c6 _PyVectorcall_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:257
19 0x00000000002173c6 _PyVectorcall_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:259
20 0x00000000002173c6 _PyObject_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:328
21 0x00000000002173c6 PyObject_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:355
22 0x00000000001eb58b do_call_core()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:7349
23 0x00000000001eb58b _PyEval_EvalFrameDefault()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:5376
24 0x000000000020d055 _PyEval_EvalFrame()  /usr/local/src/conda/python-3.11.14/Include/internal/pycore_ceval.h:73
25 0x000000000020d055 _PyEvalFrameClearAndPop()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:6406
26 0x000000000020d055 _PyEval_Vector()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:6439
27 0x000000000020d055 _PyFunction_Vectorcall()  /usr/local/src/conda/python-3.11.14/Objects/call.c:393
28 0x00000000002173c6 _PyVectorcall_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:257
29 0x00000000002173c6 _PyVectorcall_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:259
30 0x00000000002173c6 _PyObject_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:328
31 0x00000000002173c6 PyObject_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:355
32 0x00000000001eb58b do_call_core()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:7349
33 0x00000000001eb58b _PyEval_EvalFrameDefault()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:5376
34 0x000000000020d055 _PyEval_EvalFrame()  /usr/local/src/conda/python-3.11.14/Include/internal/pycore_ceval.h:73
35 0x000000000020d055 _PyEvalFrameClearAndPop()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:6406
36 0x000000000020d055 _PyEval_Vector()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:6439
37 0x000000000020d055 _PyFunction_Vectorcall()  /usr/local/src/conda/python-3.11.14/Objects/call.c:393
38 0x00000000002173c6 _PyVectorcall_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:257
39 0x00000000002173c6 _PyVectorcall_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:259
40 0x00000000002173c6 _PyObject_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:328
41 0x00000000002173c6 PyObject_Call()  /usr/local/src/conda/python-3.11.14/Objects/call.c:355
42 0x00000000001eb58b do_call_core()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:7349
43 0x00000000001eb58b _PyEval_EvalFrameDefault()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:5376
44 0x00000000002a3035 _PyEval_EvalFrame()  /usr/local/src/conda/python-3.11.14/Include/internal/pycore_ceval.h:73
45 0x00000000002a3035 _PyEval_Vector()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:6434
46 0x00000000002a277d PyEval_EvalCode()  /usr/local/src/conda/python-3.11.14/Python/ceval.c:1148
47 0x00000000002c020a run_eval_code_obj()  /usr/local/src/conda/python-3.11.14/Python/pythonrun.c:1741
48 0x00000000002bbf03 run_mod()  /usr/local/src/conda/python-3.11.14/Python/pythonrun.c:1762
49 0x00000000002bbf03 run_mod()  /usr/local/src/conda/python-3.11.14/Python/pythonrun.c:1763
50 0x00000000002d17c0 pyrun_file()  /usr/local/src/conda/python-3.11.14/Python/pythonrun.c:1657
51 0x00000000002d114a _PyRun_SimpleFileObject()  /usr/local/src/conda/python-3.11.14/Python/pythonrun.c:440
52 0x00000000002d0f24 _PyRun_AnyFileObject()  /usr/local/src/conda/python-3.11.14/Python/pythonrun.c:79
53 0x00000000002cb261 pymain_run_file_obj()  /usr/local/src/conda/python-3.11.14/Modules/main.c:360
54 0x00000000002cb261 pymain_run_file()  /usr/local/src/conda/python-3.11.14/Modules/main.c:379
55 0x00000000002cb261 pymain_run_python()  /usr/local/src/conda/python-3.11.14/Modules/main.c:605
56 0x00000000002cb261 Py_RunMain()  /usr/local/src/conda/python-3.11.14/Modules/main.c:684
57 0x0000000000292667 Py_BytesMain()  /usr/local/src/conda/python-3.11.14/Modules/main.c:738
58 0x00000000000295d0 __libc_start_call_main()  ???:0
59 0x0000000000029680 __libc_start_main_alias_2()  :0
=================================
/var/lib/slurm/slurmd/job8144292/slurm_script: line 27: 3934722 Floating point exceptionpython 456.py

crawdad · May 6, 2026, 11:24pm

I’m not yet certain what the error is, but He is a two-electron system, so there are no triple-excitation components. CCSD = FCI in this case. Also, the job appears to be triggering the numerical-differentiation driver rather than the analytic gradient code, so you wouldn’t get the density from it anyway. I don’t know how the gradient code will behave for an atom, but perhaps try CCSD (which is equivalent to FCI for He) as a test?

charlotte · May 7, 2026, 1:49am

@crawdad Thanks for the reply. Yes, CCSD works for He, and the code above also works for other atoms. However, when I try to obtain the density matrix using wfn.Da(), it returns a Nonetype object. Do you know how to get the 1RDM in psi4 instead?