Computing RDMs Segfault for FCI Class

[lcyyork]

The following input should reproduce the error.

molecule {
2 5
Fe
symmetry c1
}

set globals{
  basis                  cc-pvdz
  reference              rohf
  scf_type               df
  e_convergence          10
  d_convergence          5
  ints_tolerance         0.0
}
Escf, wfn = energy('scf', return_wfn=True)

set forte{
  int_type               df
  active_space_solver    fci
  correlation_solver     sa-mrdsrg
  frozen_docc            [5]
  restricted_docc        [4]
  active                 [5]
  e_convergence          12
  r_convergence          8
}
energy('forte', ref_wfn=wfn)

The active space solver would create three FCI objects, one for each spin Ms that is greater than 0. We will compute the RDMs for each Ms value and average them, which generates a segfault. The issue is that for Ms = 0, the max size of string list is 10; but for Ms = 2, the max size becomes 5. The energy computation starts from Ms = 0 to Ms = 1 to Ms = 2. When computing the RDMs in fci_vector_rdm.cc, a temporary matrix C is set of C1, for example, line 162. However, C1 is marked as static in line 168 of fci_vector.h, which is a 5 x 5 matrix after the Ms = 2 energy computation. Now if I want to compute the RDMs for Ms = 0, here goes the segfault.

@fevangelista So I wonder if we should make the change to FCI class or not. For my purpose (computing Ms-averaged RDMs), I can leave FCI as is and implement S+ or S- operator for the spin-complete active methods.

[fevangelista]

This should definitively be fixed. Since you have identified the main issue, it should be easy for me to solve it.

[lcyyork]

I tried to simply modify the dimension, but it got stuck for Ms = 2. I did not go further from there.

[rwashi690]

Yeah, I ran into this issue and as well with MnH septet state - here is my input and output files

import forte
molecule dimer{
0  7
Mn
H  1 R
}
dimer.R = 1.7309
set globals{
basis               cc-pVQZ
reference           rohf
scf_type            df
DOCC                [6, 0, 2, 2]
SOCC                [3, 1, 1, 1]
maxiter             250
d_convergence       1e-2
e_convergence       1e-2
damping_percentage  20
}
Escf, scf_wfn = energy('scf', return_wfn=True)
set forte {
job_type            none
avas                true
avas_diagonalize    true
subspace            ["Mn(3d)","Mn(4s)", "Mn(4p)", "H(1s)"]
avas_num_active     10
}
Enull, AVAS_wfn = energy('forte', ref_wfn=scf_wfn, return_wfn=True)
set forte{
avas                   false
subspace               []
job_type               mcscf_two_step
int_type               df
casscf_ci_solver       fci
casscf_maxiter         1000
casscf_e_convergence   1.0e-8
casscf_g_convergence   1.0e-8
casscf_micro_maxiter   75
fci_maxiter            1000
frozen_docc            [0, 0, 0, 0]
restricted_docc        [5, 0, 2, 2] # FV
active                 [5, 1, 2, 2] # FV
root_sym               0
}
Ecas, cas_wfn= energy('forte', ref_wfn=AVAS_wfn, return_wfn=True)
set forte{
job_type               newdriver
correlation_solver     sa-mrdsrg       # spin-adapted DSRG computation
corr_level             PT2          # spin-adapted theories: PT2, PT3, LDSRG2_QC, LDSRG2
avas                   false
active_space_solver    fci
subspace               []
dsrg_s                 1.5
DSRG_MAXITER           100
e_convergence          8
r_convergence          6
relax_ref              once
nroot                  1               # Only one solution
frozen_docc            [0, 0, 0, 0]
restricted_docc        [5, 0, 2, 2]    # FV : Full valence active and restricted
active                 [5, 1, 2, 2] 
root_sym               0
maxiter                1000
}
Edsrg, wfn= energy('forte', ref_wfn=cas_wfn, return_wfn=True)
rvals = [4.5, 4.4, 4.300000000000001, 4.2, 4.1, 4.0, 3.9000000000000004, 3.8000000000000003, 3.7, 3.6, 3.5, 3.4000000000000004, 3.3000000000000003, 3.2, 3.1, 3.0, 2.9000000000000004, 2.8, 2.7, 2.6, 2.5, 2.4000000000000004, 2.3, 2.2, 2.1, 2.0, 1.9, 1.8, 1.7835129715000002, 1.7624677829000002, 1.7414225943000001, 1.7309, 1.7203774057, 1.7098548114, 1.7000000000000002, 1.6993322171, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1, 1.0]
energies = {}
energies = {1.7309: Edsrg}
for r in rvals:
    dimer.R = r 
    Edsrg, wfn = energy('forte', ref_wfn=wfn, return_wfn=True)
    energies[r] = Edsrg
    psi4.print_out("\n")
    psi4.print_out("                   DSRG Final Energy                \n")
    psi4.print_out("-----------------------------------------------------\n")
    e= energies[r]
    psi4.print_out("rval: %3.10f            %10.20f\n" % (r, e))

Output_Fail.TXT

[rwashi690]

Also MnH Quintet runs perfectly fine

[fevangelista]

I posted the fix in a new PR (#269). @rwashi690: once merged, can you make sure that you can run MnH septet?