Help with FCIMQC Setup for Cyclo[16]carbon Using NECI — Different Results from SCF vs CASSCF Orbitals

[Marcovit94]

Hi there,

I'm a beginner working with QMC and trying to set up a FCIMQC calculation using NECI for a C₁₆ cyclocarbon, focusing on the full π-electron active space (32 electrons in 16 orbitals). I've generated two different FCIDUMP files from OpenMolcas and get quite different results when using them in NECI.

Approach 1: SCF orbitals

&GATEWAY
 COORD = opted.xyz
 title = test
 group = z
 basis = cc-pvdz

&SEWARD
&SCF
&RASSCF
 Spin = 1
 CIRoot = 1 1 1
 Frozen = 16 0
 Inactive = 16 0
 RAS1 = 0 0
 RAS2 = 16 16
 RAS3 = 0 0
 NActEl = 32 0 0
 DMPO

This setup produces a cleanly converging NECI run with a total energy of: E = -605.72129172258 Ha

Approach 2: CASSCF orbitals

&GATEWAY
 COORD = opted.xyz
 title = test
 group = z
 basis = cc-pvdz

&SEWARD
&SCF
&RASSCF
 Spin = 1
 CIRoot = 1 1 1
 Frozen = 16 0
 Inactive = 21 5
 RAS1 = 0 0
 RAS2 = 6 6
 RAS3 = 0 0
 NActEl = 12 0 0

&RASSCF
 Spin = 1
 CIRoot = 1 1 1
 Frozen = 16 0
 Inactive = 16 0
 RAS1 = 0 0
 RAS2 = 16 16
 RAS3 = 0 0
 NActEl = 32 0 0
 DMPO

This was an attempt to improve the orbital quality by first doing a smaller π-π* CASSCF (12,12) and then using those orbitals for the full space. However, with this FCIDUMP, NECI gives worse and unstable results. The best I got so far is: E = -604.54366646533 Ha
This is more than 1.2 Ha higher than the previous result, and I'm not sure why. I have also tried (4,4) and (8,8) active spaces and got the same result.

NECI Input

Here's the NECI input I use for both:

system read
    nonuniformrandexcits pchb delocalised
    electrons 32
endsys

calc
    methods
        method vertex fcimc
    endmethods

    load-balance-interval 1000
    nmcyc              5000
    totalwalkers       10000000

    diagshift           0.0
    startsinglepart      10
    truncinitiator
    addtoinitiator        3
    allrealcoeff
    no-changeref
    realspawncutoff     0.7
    avgrowthrate
    shiftdamp           0.1
    stepsshift            1
    memoryfacpart       5.0
    memoryfacspawn     10.0
    tau-values \
        start user-defined 0.02 \
        max 0.05
    tau-search \
        algorithm conventional \
        stop-condition n-opts 5 \
        maxwalkerbloom 1
endcalc

end

i also tried nonuniformrandexcits pchb localised but have the same exact issue.

My question(s):

1. Why might the CASSCF-generated orbitals yield significantly worse FCIMC performance or higher energy than plain SCF orbitals?
2. Is there something I should change in the NECI input to better accommodate casscf orbitals or improve convergence?
3. Is it common that orbitals from a smaller CASSCF don't help, or possibly hurt, for larger full-valence FCI-type runs?

I’d really appreciate any insights — especially as I’m just starting out with QMC and NECI. Thanks in advance!

I've attached the FCIMCStats output files from both NECI runs for comparison.

FCIMCStats_SCF.txt

FCIMCStats_CASSCF.txt

Best,
Marco

[MaruSongMaru]

Hi Marco, thank you for posting the issue.

Based on the NECI input and stat files you shared, I speculate that the reference determinant chosen by NECI of the CASSCF calculation (or maybe both calculations) is not the correct one. Ideally, the reference determinant has to be the leading determinant in the exact wave function. In delocalized basis, the Hartree-Fock determinant is the one in general.

In the stat files, I can see the projected energy (column 23) at iteration 0 of CASSCF is greater than SCF, meaning that the energy of the reference determinant of CASSCF is greater than SCF. Also, reference determinant population (column 12) of CASSCF is decreasing to zero, suggesting that the reference determinant of the calculation may not be the actual leading determinant. As a result, you can see the projected energy of CASSF is very noisy, as the denominator of projected energy is the reference population.

To fix this issue, I would suggest to modify the NECI input file as follows:

1. Explicitly set the reference determinant using definedet keyword in the calc block. A set of integers follow definedet in the input. The integers are spin-orbital indices (odd indices are beta electrons and even indices are alpha electrons). For example, if orbitals in FCIDUMP are sorted by energy and you run a 4-electron calculation, then definedet 1 2 3 4 sets the HF determinant (1/2 are beta/alpha electrons in orbital 1, and so on). In your case, definedet 1-32 might be the HF if the orbitals are sorted by energy.
   For details of the keyword, please see the documentation.
2. Remove no-changeref keyword. NECI, by default, updates reference determinant based on instantaneous walker populations on determinants. With no-changeref, the reference determinant is fixed. You can use proje-changeRef keyword to control the sensitivity of reference update.
3. (optional) I recommend to use semistochastic feature to reduce stochastic noise. You can add, for example,

semi-stochastic 2000
pops-core 100000

in the calc block.

To summarize, I guess it is a problem of a bad reference choice. If that is the case, adding the following keywords in the calc block can help.

calc
...
definedet 1-32 # please change the numbers according to your FCIDUMP
# no-changeref

semi-stochastic 2000
pops-core 100000
...
endcalc

[glimanni]

I might be wrong but I also suspect some error in the preceding cas1212. Is that calculation lower or higher in energy than your SCF? Maybe your choice of frozen/inactive is such that you converge to an undesired ion? Let us know if that is of any help. Giovanni

…

On Wed, Jun 4, 2025 at 17:55 Maru Song ***@***.***> wrote: *MaruSongMaru* left a comment (fkfest/NECI_STABLE#19) <#19 (comment)> Hi Marco, thank you for posting the issue. Based on the NECI input and stat files you shared, I speculate that the reference determinant chosen by NECI of the CASSCF calculation (or maybe both calculations) is not the correct one. Ideally, the reference determinant has to be the leading determinant in the exact wave function. In delocalized basis, the Hartree-Fock determinant is the one in general. In the stat files, I can see the projected energy (column 23) at iteration 0 of CASSCF is greater than SCF, meaning that the energy of the reference determinant of CASSCF is greater than SCF. Also, reference determinant population (column 12) of CASSCF is decreasing to zero, suggesting that the reference determinant of the calculation may not be the actual leading determinant. As a result, you can see the projected energy of CASSF is very noisy, as the denominator of projected energy is the reference population. To fix this issue, I would suggest to modify the NECI input file as follows: 1. Explicitly set the reference determinant using definedet keyword in the calc block. A set of integers follow definedet in the input. The integers are spin-orbital indices (odd indices are beta electrons and even indices are alpha electrons). For example, if orbitals in FCIDUMP are sorted by energy and you run a 4-electron calculation, then definedet 1 2 3 4 sets the HF determinant (1/2 are beta/alpha electrons in orbital 1, and so on). In your case, definedet 1-32 might be the HF if the orbitals are sorted by energy. For details of the keyword, please see the documentation <https://www2.fkf.mpg.de/alavi/neci/stable/page/01_user_doc/03_calculation_inputs.html> . 2. Remove no-changeref keyword. NECI, by default, updates reference determinant based on instantaneous walker populations on determinants. With no-changeref, the reference determinant is fixed. You can use proje-changeRef keyword to control the sensitivity of reference update. 3. (optional) I recommend to use semistochastic feature to reduce stochastic noise. You can add, for example, semi-stochastic 2000 pops-core 100000 in the calc block. To summarize, I guess it is a problem of a bad reference choice. If that is the case, adding the following keywords in the calc block can help. calc ... definedet 1-32 # please change the numbers according to your FCIDUMP # no-changeref semi-stochastic 2000 pops-core 100000 ... endcalc — Reply to this email directly, view it on GitHub <#19 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB6APTXQRUQSFCWYWMIRRRL3B4JFZAVCNFSM6AAAAAB6SUWS72VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDSNBQGU3DENRUGA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[Marcovit94]

Hi all – thanks again for the very detailed advice; it has already helped a lot.
Below are the benchmarks I ran after applying your suggestions, plus a couple of follow-up questions about definedet and restarts.

CASSCF(12,12) sanity check:
I followed this paper on cyclo-[16]carbon which recommends a balanced 12,12 in-plane vs out-of-plane π active space: https://pubs.acs.org/doi/10.1021/jacs.3c10207

Energies (Molcas, cc-pVDZ, same geometry)
RHF              –605.324460270
CASSCF(12,12)    –605.40838957

So the CASSCF is ~0.084 Ha lower than RHF, which looks reasonable and gives me confidence the active space isn’t pathological—but please correct me if I’m mis-interpreting, I’m still learning all this.

FCIQMC results with your suggested settings:
I re-ran NECI with:

# removed no-changeref
semi-stochastic 2000
pops-core       100000

and progressively increased the number of walkers (1 M → 5 M). The projected energy stabilises nicely; for example, one of the 5 M walker runs gives: Final projected energy ≈ –605.76055644347 Ha

So this looks much better!

Confusion about definedet:
I tried setting definedet 1-32 manually, but the reference determinant started at: E = –600.312280532 Ha which is ~5 Ha above RHF, so that’s clearly a bad reference. The run doesn't converge well from this. Then I tried not setting definedet, and just let NECI find the best reference. That time, it eventually switched to: 1-16 33-48 which had energy of E = –605.316496600 Ha which is much more reasonable, close to RHF energy but its still higher than the casscf energy CASSCF(12,12) –605.40838957. That run converges smoothly and ends up around the qmc energy: ≈ –605.760

So clearly my manual 1-32 was wrong for this FCIDUMP, but I’m still not sure how to choose the correct definedet manually from the beginning.

I’m guessing it’s related to this part of the Molcas log?

printout of CI-coefficients larger than 0.05 for root 1
energy= -605.408390
conf/sym  111111 222222     Coeff  Weight
   22718  222000 222000   0.93144 0.86759
   22722  222000 220200  -0.05275 0.00278

The leading 222000 222000 occupation obviously isn’t just 16 × 2, so I’m not sure how to map that into the correct spin-orbital indices for definedet. Any explanation on how to decode this (or a general rule) would be super helpful.

Restart question:
When I restart a QMC run, is it enough to just use readpops, or should I also reset definedet to match the last one that was being used (e.g. 1-16 33-48 in this case)? I’m not sure if it matters once walkers are loaded.

Files attached:
I’ve attached the following:

molcas.log

FCIDUMP

opted.xyz (geometry)

NECI logs:

    good run (auto reference): neci_good.out, 

    bad run (definedet 1-32): neci_bad.out, FCIMCStats_bad

I’d really like to understand how to set up definedet properly, especially because I want to scan through different ring sizes of cyclocarbons next. Thanks again for all the help and detailed answers — it’s super appreciated!

Best,
Marco

molcas.log.txt

opted.xyz.txt

FCIDUMP.txt

FCIMCStats_bad.txt

neci_bad.out.txt

neci_good.out.txt

[glimanni]

You are correct about definedet setup. The key is in Molcas: 222000,222000. Molcas uses space orbitals, while NECI uses Spin-orbitals. So 2 spin orbitals map one space orbitals. Moreover, you are using point group symmetry. So you have some doubly occupied and some empty in each Irrep. The correct definedet should be: 1,2,3,4,5,6, 13,14,15,16,17,18 Spin orbitals from 7-12 are empty so don’t show up in definedet. In the definedet you may read odd numbers as spin-orbital with beta electrons and even numbers as spin-orbs with Alpha electrons. Hope that helps you on that part. Giovanni

…

On Fri, Jun 6, 2025 at 12:22 Marco Vitek ***@***.***> wrote: *Marcovit94* left a comment (fkfest/NECI_STABLE#19) <#19 (comment)> Hi all – thanks again for the very detailed advice; it has already helped a lot. Below are the benchmarks I ran after applying your suggestions, plus a couple of follow-up questions about definedet and restarts. *CASSCF(12,12) sanity check:* I followed this paper on cyclo-[16]carbon which recommends a balanced 12,12 in-plane vs out-of-plane π active space: https://pubs.acs.org/doi/10.1021/jacs.3c10207 Energies (Molcas, cc-pVDZ, same geometry) RHF –605.324460270 CASSCF(12,12) –605.40838957 So the CASSCF is ~0.084 Ha lower than RHF, which looks reasonable and gives me confidence the active space isn’t pathological—but please correct me if I’m mis-interpreting, I’m still learning all this. *FCIQMC results with your suggested settings:* I re-ran NECI with: # removed no-changeref semi-stochastic 2000 pops-core 100000 and progressively increased the number of walkers (1 M → 5 M). The projected energy stabilises nicely; for example, one of the 5 M walker runs gives: Final projected energy ≈ –605.76055644347 Ha So this looks much better! *Confusion about definedet:* I tried setting definedet 1-32 manually, but the reference determinant started at: E = –600.312280532 Ha which is ~5 Ha above RHF, so that’s clearly a bad reference. The run doesn't converge well from this. Then I tried not setting definedet, and just let NECI find the best reference. That time, it eventually switched to: 1-16 33-48 which had energy of E = –605.316496600 Ha which is much more reasonable, close to RHF energy but its still higher than the casscf energy CASSCF(12,12) –605.40838957. That run converges smoothly and ends up around the qmc energy: ≈ –605.760 So clearly my manual 1-32 was wrong for this FCIDUMP, but I’m still not sure how to choose the correct definedet manually from the beginning. I’m guessing it’s related to this part of the Molcas log? printout of CI-coefficients larger than 0.05 for root 1 energy= -605.408390 conf/sym 111111 222222 Coeff Weight 22718 222000 222000 0.93144 0.86759 22722 222000 220200 -0.05275 0.00278 The leading 222000 222000 occupation obviously isn’t just 16 × 2, so I’m not sure how to map that into the correct spin-orbital indices for definedet. Any explanation on how to decode this (or a general rule) would be super helpful. *Restart question:* When I restart a QMC run, is it enough to just use readpops, or should I also reset definedet to match the last one that was being used (e.g. 1-16 33-48 in this case)? I’m not sure if it matters once walkers are loaded. *Files attached:* I’ve attached the following: molcas.log FCIDUMP opted.xyz (geometry) NECI logs: good run (auto reference): neci_good.out, bad run (definedet 1-32): neci_bad.out, FCIMCStats_bad I’d really like to understand how to set up definedet properly, especially because I want to scan through different ring sizes of cyclocarbons next. Thanks again for all the help and detailed answers — it’s super appreciated! Best, Marco molcas.log.txt <https://github.com/user-attachments/files/20626252/molcas.log.txt> opted.xyz.txt <https://github.com/user-attachments/files/20626255/opted.xyz.txt> FCIDUMP.txt <https://github.com/user-attachments/files/20626261/FCIDUMP.txt> FCIMCStats_bad.txt <https://github.com/user-attachments/files/20626262/FCIMCStats_bad.txt> neci_bad.out.txt <https://github.com/user-attachments/files/20626264/neci_bad.out.txt> neci_good.out.txt <https://github.com/user-attachments/files/20626269/neci_good.out.txt> — Reply to this email directly, view it on GitHub <#19 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB6APTVE4L2CG7M3VHR45AL3CFTWRAVCNFSM6AAAAAB6SUWS72VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDSNBYHAYDCNZYGM> . You are receiving this because you commented.Message ID: ***@***.***>

[MaruSongMaru]

definedet choice
Choosing a good reference determinant of course depends on how orbitals in the FCIDUMP file are ordered/labeled.
Since I'm not an expert on Molcas, I just guess based on your results that the first 16 orbitals in FCIDUMP belong to the sym 1 sector and the next 16 orbitals belong to the sym 2 sector and they are individually sorted by orbital energy, so I guess that is why 1-16 33-48 (spatial orbital 1 - 8 and 17 - 24) is a good reference.

The reason why your "bad" calculation is bad is probably because the total population is only 1e5, not because of the bad definedet of the bad calculation.
As you can see in the output of the bad calculation, NECI dynamically updates the reference determinant (grep Ref) based on the instantaneous walker distribution. In the later part of the output, you can also see the instantaneous walker distribution at the last iteration of the calculation:

 Excitation   ExcitLevel   Seniority    Walkers    Amplitude    Init?   <D|H|D>  Proc  
1111111111111110111010000100000011110111110111000000000000000000    7    6   -31.68900      0.09608  Y  -6.04075559E+02      31    
1111111111111111111010000000000011111111110101000000000000000000    7    4   -23.92661      0.07254  Y  -6.04582591E+02      15    
1111111111110111111011000000000011111111110011000000000000000000    6    2   -19.80419      0.06004  Y  -6.04230130E+02       8     
1111111111111111111001000000000011111111110011000000000000000000    6    2   -11.63276      0.03527  Y  -6.04546973E+02      19    

Each line is a determinant and they are sorted by the instantaneous walker distribution. The first column is a bitstring representation of spin orbitals. 1 means occupied and 0 means unoccupied. As you can see, the wavefunction is quite multi reference (the third column is the number of walkers on each determinant). If you see this part of the good calculation output, you can find a relatively compact wave function. This suggests that your total population of the calculation (5M) is enough to solve the problem, while 1e5 is not enough.
This also means, if you set the target population enough to solve the problem, the correct wave function "emerges" during the course of the simulation and the reference determinant is updated accordingly, regardless of the initial reference.

Restart
As far as I know, a popsfile does not contain the reference determinant information. So, the reference choice at the initialization stage is not affected by the popsfile. While initializing the calculation, NECI sets a reference determinant either using its own algorithm or using definedet (if provided).
But, regardless of the initial reference, NECI immediately updates the reference based on the initial walker distribution given by the popsfile, unless you have no-changeref.

NECI tip
In case you just use the final energy written in the NECI output, I would recommend not to use the energy because iterations before the dynamics is stabilized can be included.
Instead, I would suggest to do error analysis on the stat file using the iterations after convergence.
There is a python script for blocking analysis in the NECI repo (see https://www2.fkf.mpg.de/alavi/neci/stable/page/01_user_doc/08_error_analysis.html).

[Marcovit94]

Thank you both for all the help—this thread has been super useful, and I think I understand what’s going on with definedet, walker populations, and restarts.

A more general question on scalability

I’m now looking ahead to larger systems and haven’t found a clear answer in the literature:

I have two π-isomers, each with 88 π-electrons in nearly identical orbital manifolds. How plausible is it to treat the full π system—i.e. an (88,88) active space—with NECI without monopolising an entire cluster for weeks?

My thinking is:

• Absolute energies at this size likely won’t converge without an enormous number of walkers.

• But maybe relative energies between the isomers could converge faster due to error cancellation?

Or is that still too optimistic, and I should instead define a smaller active space?

The largest active space I found in the literature was CAS(32,34)—so I’m curious if anyone has pushed NECI beyond that, and whether it worked for relative comparisons.

Any input would be really appreciated!

Thanks again,
Marco

[glimanni]

Good observations! It is true that you won’t be able to converge the total energies. I’d strongly suggest to switch to stochastic-GAS as discussed in here: Stochastic Generalized Active Space Self-Consistent Field: Theory and Application | Journal of Chemical Theory and Computation https://pubs.acs.org/doi/10.1021/acs.jctc.1c00936 It is relatively recent and we used it for cases similar to the one you are describing. There you can get inspiration for a valuable setups for your system. Inputs are available in the SI. If you need more information please do not hesitate to contact us. If you would like to discuss your options we could setup a zoom meeting together. 😉 Best, Giovanni

…

On Thu, Jun 19, 2025 at 15:48 Marco Vitek ***@***.***> wrote: *Marcovit94* left a comment (fkfest/NECI_STABLE#19) <#19 (comment)> Thank you both for all the help—this thread has been super useful, and I think I understand what’s going on with definedet, walker populations, and restarts. *A more general question on scalability* I’m now looking ahead to larger systems and haven’t found a clear answer in the literature: I have two π-isomers, each with 88 π-electrons in nearly identical orbital manifolds. How plausible is it to treat the full π system—i.e. an (88,88) active space—with NECI without monopolising an entire cluster for weeks? My thinking is: - Absolute energies at this size likely won’t converge without an enormous number of walkers. - But maybe relative energies between the isomers could converge faster due to error cancellation? Or is that still too optimistic, and I should instead define a smaller active space? The largest active space I found in the literature was CAS(32,34)—so I’m curious if anyone has pushed NECI beyond that, and whether it worked for relative comparisons. Any input would be really appreciated! Thanks again, Marco — Reply to this email directly, view it on GitHub <#19 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB6APTWYQTXE7XSBDY34ON33EK5TZAVCNFSM6AAAAAB6SUWS72VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDSOBYGE3TIMBSGQ> . You are receiving this because you commented.Message ID: ***@***.***>