Open shell calculation of Li. Strange RDM1 occupation with half of electron and n_alpha - n_beta=0

[dithillobothrium]

Using the FCIDUMP generated by the method proposed from previous issue I'm trying to calculate Li atom with 3 electrons.

Input file is

num_thrds 12
nelec 3
orbitals li_6orb.fcidump
spin   1

schedule default
maxiter 20
startM 100
maxM 1000
hf_occ 2 1 0 0 0 0

onepdm

output reports that Spin=1

Finished Sweep with 1000 states and sweep energy for State [ 0 ] with Spin [ 1 ] :: -7.395656723881

but density matrix onepdm.0.0.txt has following occupations

12
0 0 9.99984987454079e-01
1 1 9.99984987454079e-01
2 2 4.99992968823474e-01
2 3 0.00000000000000e+00
3 2 0.00000000000000e+00
3 3 4.99992968823474e-01

So, the non-paired electron is divided between spin-up and spin-down orbital and non-diagonal components are zero. The number of electrons for spin-up is 1.5, and n_alpha - n_beta = 0.

Why does it happen? I'm not sure it works correctly.

Best regards

[dithillobothrium]

Sorry, I need to assign someone here @sanshar @gkc1000 @hczhai

[hczhai]

Since you do not have nonspinadapted in the input file, this is SU(2) spin-adapted calculation (alpha = beta).
In the input file spin 1 means total spin 2S = 1.
If you want alpha and beta to be different, please add nonspinadapted in the input file.

[gkc1000]

For computational efficiency the code uses the singlet embedding strategy to treat nonsinglet states with SU2 as described in Sharma and Chan. Additional fictitious spins are added so that the final state is a singlet. This means that the ensemble average over all Ms sectors is represented rather than a given Ms sector. I think it is possible to turn this off as an input option but I think the two particle DM might not have been coded correctly without it. Sandeep is this true? Huanchen can you check the input options?

On Thu, Feb 18, 2021 at 5:28 AM hczhai ***@***.***> wrote: Since you do not have nonspinadapted in the input file, this is SU(2) spin-adapted calculation (alpha = beta). In the input file spin 1 means total spin 2S = 1. If you want alpha and beta to be different, please add nonspinadapted in the input file. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#56 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN5NRAY3ZOWVQWYWFZLIWLS7UIWDANCNFSM4XLGKXFA> .

-- Sent from Gmail Mobile

[dithillobothrium]

Since you do not have nonspinadapted in the input file, this is SU(2) spin-adapted calculation (alpha = beta).
In the input file spin 1 means total spin 2S = 1.
If you want alpha and beta to be different, please add nonspinadapted in the input file.

@hczhai
for a total spin it should be 2S+1 = 2 in my case. From the manual I understood I should put a difference in electrons
as na-nb=1. Spin adaptation means making a state which is eigenstate of S^2. But, as I remember Sz should give a na-nb=1 in any case. How can I reuse this 1RDM after the calculation? I would like to upload it back in the DFT/HF calculation and plot density. What is a basis of RDM ? How can I use the molecular orbitals/coefficients for that, which were used for generation of FCIDUMP?

For computational efficiency the code uses...

@gkc1000
How physically correct is this approach? The same question, how can I reuse this RDM? Should I transform it somehow? Could you please send a link to this paper you mentioned?

Thank you.

[dithillobothrium]

You mean the final rdm is rdm taken from a mixed state?
In this case rdm becomes not very useful, since it is averaged over ensemble. But I don't understand connection between adding fictious spin and making a mixed state.

For computational efficiency the code uses the singlet embedding strategy to treat nonsinglet states with SU2 as described in Sharma and Chan. Additional fictitious spins are added so that the final state is a singlet. This means that the ensemble average over all Ms sectors is represented rather than a given Ms sector. I think it is possible to turn this off as an input option but I think the two particle DM might not have been coded correctly without it. Sandeep is this true? Huanchen can you check the input options?
On Thu, Feb 18, 2021 at 5:28 AM hczhai @.***> wrote: Since you do not have nonspinadapted in the input file, this is SU(2) spin-adapted calculation (alpha = beta). In the input file spin 1 means total spin 2S = 1. If you want alpha and beta to be different, please add nonspinadapted in the input file. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#56 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRAY3ZOWVQWYWFZLIWLS7UIWDANCNFSM4XLGKXFA .
-- Sent from Gmail Mobile

[hczhai]

@gkc1000

I think it is possible to turn this off as an input option but I think the two particle DM might not have been coded correctly without it. Sandeep is this true? Huanchen can you check the input options?

(a) The singlet embedding is controlled by the internal variable m_add_noninteracting_orbs = true (https://github.com/sanshar/StackBlock/blob/master/input.C#L140). Unfortunately, this parameter is not linked with any parameter in the input file, in the most recent version of StackBlock.
(b) In a fork of the old Block for spin-orbit coupling (https://github.com/ElviraRS/Block), there is an input parameter to disable singlet embedding: https://github.com/ElviraRS/Block/blob/master/input.C#L581-L586 However, without singlet embedding, the 1pdm will still be same, namely, for the mixed state.
(c) In the StackBlock code, for the spin-adapted mode, the alpha/beta components of onepdm is explicitly coded as half of the "spatial" onepdm: https://github.com/sanshar/StackBlock/blob/master/modules/onepdm/onepdm.C#L136-L143. Therefore, it will always be alpha = beta.
(d) The correct 1pdm for M != 0 can be obtained from nonspinadapted mode.

[hczhai]

@dithillobothrium

for a total spin it should be 2S+1 = 2 in my case. From the manual I understood I should put a difference in electrons
as na-nb=1. Spin adaptation means making a state which is eigenstate of S^2.

In StackBlock, when you are in spin-adapted mode, the input parameter spin has the meaning of "total spin", which can only be a non-negative number. When you are in non-spin-adapted mode, the input parameter spin has the meaning of "projected spin", which can be positive or negative, or zero.

But, as I remember Sz should give a na-nb=1 in any case.

It can also be nb-na=1. The code does not know which case it is. In the onepdm you get a mixture of the two cases.