Incorrect number of electrons and magnetizabilities

[ejberquist]

This might be two bugs, but one led to the other.

1. I'm trying to run magnetizability calculations, testing a simple system just to get some numbers (LiH+ doublet, PBE/3-21G), but the program crashes in ABACUS:

Code: Select all

   Number of electrons from numerical integration:      1.999999
   Number of electrons from orbital occupations:        3.000000
   Error in the number of electrons:                   -1.000001
   Error larger than DFTELS (set input):                0.001000
   Calculation aborted.

This made me think that it could be due to DFT, the small basis set, or an occupation problem. Increasing the basis set to cc-pVDZ doesn't change anything. Switching from PBE to HF causes the calculation to complete, because there isn't a check for electron count, but I'm assuming there is still an issue. To confirm this isn't an occupation problem, I tried setting the number of electrons explicitly, which didn't alter the result, but setting

Code: Select all

*SCF INPUT
.DOUBLY OCCUPIED
1
.SINGLY OCCUPIED
1

for HF changed two of the magnetizability principal components by roughly a factor of two. The PBE calculation still crashed due to the electron count check.

2. The number of electrons are reported incorrectly for any open-shell HF or DFT calculation; active electrons aren't counted, only 2*NISHT. From the same calculation, which should be a doublet with 3 electrons:

Code: Select all

 Number of electrons :    2
 Orbital occupations :    1

and for the triplet with 4 electrons:

Code: Select all

 Number of electrons :    2
 Orbital occupations :    1

The Wave function specification block is correct, though.

[taylor]

I cannot comment on open-shell DFT property calculations as I very rarely do them, but I think there is a requirement that (unlike open-shell HF) they have to be run using

Code: Select all

.DIRECT

(or .PARALLEL).

You can easily check whether your doublet HF calculation is correct: run a CASSCF calculation with a single open shell and one electron. I am sure this exercises the MCSCF code, not SCF and this has been used many times over the years and I would be very surprised if there is a bug in it. If you want to be absolutely sure it is MCSCF, create a three active orbital space from the sigma and the two pi orbitals with the one active electron. Symmetry will guarantee that there is only one configuration possible, although you might be wise in that case to start using symmetry explicitly so you can see check is going on...

Best regards
Pete

[taylor]

Garbled the last edit which was to say "so you can check to see what is going on". It may for instance be that running without symmetry and selecting three active orbitals will not be one sigma and two pi's, but two sigma's and one pi, which would be broken symmetry and create other problems. Run with symmetry. You will get 4 irreps and the INACTIVE will be 1 0 0 0 and the active would be 1 1 1 0, with one active electron.

Best regards
Pete

[olav]

To my knowledge, properties calculated with London orbitals have not been implemented for open-shell HF/DFT. That is, you need to use the MCSCF module to calculate magnetizabilities of open-shell systems.

Regards,
Olav

[ejberquist]

Thank you everyone. There are certainly no issues with the energies or occupations, as the ROHF calculations with and without symmetry give the same result as MCSCF active spaces of [1,0,0,0] and [1,1,1,0] for 1 electron and [2,0,0,0] for 3 electrons, all with symmetry, as expected. I don't think I'd ever want to run an MCSCF calculation *without* symmetry on a symmetric system...

Though, there are still some strange issues. Both the ROHF and [2,0,0,0]/3 give identical magnetizability tensors out to the 9th or 10th decimal place, which I think is to be expected in the sum-over-states picture, but this partially goes against Olav's statement about the London orbitals. The results are very different when not using London orbitals, but still identical between ROHF and MCSCF. Perhaps open-shell DFT isn't working for some other reason.

The other thing I don't understand is that the [1,0,0,0]/1 and [1,1,1,0]/1 calculations crash after SIRIUS, and the errors aren't clear to me.

[olav]

I don't understand, you say "ROHF calculations with and without symmetry give the same result as MCSCF active spaces of [1,0,0,0] and [1,1,1,0] for 1 electron", yet these are the attached calculations that crash.

Please attach the [2,0,0,0]/3 output as well

/Olav

[ejberquist]

I should have been more clear; they all give the same energy, but those two crash trying to calculate the magnetizability, while the [2,0,0,0]/3 calculation is successful.