qiskit-community / qiskit-nature

Qiskit Nature is an open-source, quantum computing, framework for solving quantum mechanical natural science problems.
https://qiskit-community.github.io/qiskit-nature/
Apache License 2.0
301 stars 206 forks source link

Strange behavior of $S^2$ expectation values #816

Closed MartinBeseda closed 2 years ago

MartinBeseda commented 2 years ago

Environment

What is happening?

While testing the expectation values of $S^2$ operator, it seems, that for singlet and triplet states there are expectation values 2 and 0, respectively, while the opposite way would be expected.

The states tested:

$|\psi_1 \rangle = \frac{1}{\sqrt{2}}(|0110\rangle - |1001\rangle)$ $|\psi_2 \rangle = \frac{1}{\sqrt{2}}(|0110\rangle + |1001\rangle)$

How can we reproduce the issue?

Both the cases with their expectation values are demonstrated when running the following code:

#!/usr/bin/env python3
import qiskit_nature
from qiskit.providers.aer import Aer
from qiskit.quantum_info import Statevector
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit_nature.drivers import Molecule
from qiskit_nature.drivers.second_quantization import ElectronicStructureMoleculeDriver, ElectronicStructureDriverType
from qiskit_nature.mappers.second_quantization import JordanWignerMapper
from qiskit_nature.problems.second_quantization import ElectronicStructureProblem
from qiskit_nature.transformers.second_quantization.electronic import ActiveSpaceTransformer

if __name__ == '__main__':
    qiskit_nature.settings.dict_aux_operators = True
    backend = Aer.get_backend('statevector_simulator')
    geometry = [('N', [0.000000000000, 0.000000000000, 0.000000000000]),
                ('C', [0.000000000000, 0.000000000000, 1.498047000000]),
                ('H', [0.000000000000, -0.938765985000, 2.004775984000]),
                ('H', [0.000000000000, 0.938765985000, 2.004775984000]),
                ('H', [-0.744681452, -0.131307432, -0.634501434])]

    # Obtaining ground state
    driver = ElectronicStructureMoleculeDriver(Molecule(geometry),
                                               basis="sto3g",
                                               driver_type=ElectronicStructureDriverType.PSI4)
    as_transformer = ActiveSpaceTransformer(num_electrons=2, num_molecular_orbitals=2)
    es_problem = ElectronicStructureProblem(driver, transformers=[as_transformer])
    converter = QubitConverter(JordanWignerMapper())

    # Specify statevectors
    statevec = Statevector([0, 0, 0, 0, 0, 0, 0.70710678, 0, 0, -0.70710678, 0, 0, 0, 0, 0, 0])
    statevec.draw('latex')

    statevec2 = Statevector([0, 0, 0, 0, 0, 0, 0.70710678, 0, 0, 0.70710678, 0, 0, 0, 0, 0, 0])
    statevec2.draw('latex')

    # Obtain S^2 operator
    ss_op = es_problem.second_q_ops()['AngularMomentum']
    print(f'S^2 operator:\n{ss_op}')

    ss_op_qubit = converter.convert(ss_op)
    ss_op_mat = ss_op_qubit.to_matrix().real
    print(f'S^2 matrix:\n{ss_op_mat}')

    # Obtain expectation values
    print(f'<psi|S^2|psi> = {statevec.data @ ss_op_mat @ statevec.data}')
    print(f'<psi2|S^2|psi2> = {statevec2.data @ ss_op_mat @ statevec2.data}')

Output:

S^2 operator:
Fermionic Operator
register length=4, number terms=12
  (0.75+0j) * ( +_0 -_0 )
+ (0.75+0j) * ( +_1 -_1 )
+ (0.75+0j) * ( +_2 -_2 )
+ (0.75+0j) * ( +_3 -_3 )
+ (0.5+0j) * ( +_0 -_0 +_1 -_1 )
+ (-1.5+0j) * ( +_0 -_0 +_2 -_2 )
+ (-0.5+0j) * ( +_0 -_0 + ...
S^2 matrix:
[[ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.75  0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.75  0.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    2.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.75  0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    1.    0.    0.   -1.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.75  0.    0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.75  0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.   -1.    0.    0.    1.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.75
   0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
   2.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.75  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.75  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
   0.    0.    0.    0.  ]]
<psi|S^2|psi> = (1.9999999932878738+0j)
<psi2|S^2|psi2> = 0j

What should happen?

I'd expect for the singlet state, i.e. $|\psi_1\rangle$ to be associated with the expectation value 0 and for $|\psi_2\rangle$ to be associated with 2.

Any suggestions?

No response

mrossinek commented 2 years ago

I think you are confusing $\psi_1$ and $\psi_2$ here. Let me try to show you why:

The following script, is an adaptation of yours, which uses the ExcitedStatesEigensolver in combination with a NumPyEigensolver to compute the ground and all excited states of your system. It then uses the obtained Statevector objects to evaluate all observables (the solver does this already, but I want to print things the way you did).

import numpy as np
from qiskit_nature.algorithms.excited_states_solvers import (
    ExcitedStatesEigensolver, NumPyEigensolverFactory)
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit_nature.drivers import Molecule
from qiskit_nature.drivers.second_quantization import (
    ElectronicStructureDriverType, ElectronicStructureMoleculeDriver)
from qiskit_nature.mappers.second_quantization import JordanWignerMapper
from qiskit_nature.problems.second_quantization import \
    ElectronicStructureProblem
from qiskit_nature.transformers.second_quantization.electronic import \
    ActiveSpaceTransformer
from qiskit_nature.settings import settings

if __name__ == "__main__":
    settings.dict_aux_operators = True
    geometry = [
        ("N", [0.000000000000, 0.000000000000, 0.000000000000]),
        ("C", [0.000000000000, 0.000000000000, 1.498047000000]),
        ("H", [0.000000000000, -0.938765985000, 2.004775984000]),
        ("H", [0.000000000000, 0.938765985000, 2.004775984000]),
        ("H", [-0.744681452, -0.131307432, -0.634501434]),
    ]

    # Obtaining ground state
    driver = ElectronicStructureMoleculeDriver(
        Molecule(geometry), basis="sto3g", driver_type=ElectronicStructureDriverType.PYSCF
    )
    as_transformer = ActiveSpaceTransformer(num_electrons=2, num_molecular_orbitals=2)
    es_problem = ElectronicStructureProblem(driver, transformers=[as_transformer])
    converter = QubitConverter(JordanWignerMapper())

    q_ops = {name: converter.convert(op) for name, op in es_problem.second_q_ops().items()}

    algo = ExcitedStatesEigensolver(converter, NumPyEigensolverFactory())
    res = algo.solve(es_problem)
    for val, vec in zip(res.eigenenergies, res.eigenstates):
        state = vec.primitive
        cleaned_state_dict = {k: v for k, v in state.to_dict().items() if not np.abs(v) < 1e-6}
        print(f"State: {cleaned_state_dict}")
        for name, q_op in q_ops.items():
            mat = q_op.to_matrix().real

            res = state.data @ mat @ state.data
            print(f"\t<psi|{name}|psi> = {0.0 if np.abs(res.real) < 1e-6 else res.real}")
        print()
Expand this for the full output ``` State: {'0101': (0.943726286909139+0j), '0110': (-0.169197999288364+0j), '1001': (-0.16919799928804474+0j), ' 1010': (-0.22830849627297578+0j)} = 2.000000000000001 = -1.1709917394922589 = 0.0026004017257202448 = 0.02908727852216155 = 1.3957696360113352 = 0.0 = 0.0 State: {'0011': (1+0j)} = 2.0 = -1.169966734319283 = 0.13275615748360492 = -0.04702411665523931 = 2.0316863985002866 = 2.0 = 1.0 State: {'1100': (1+0j)} = 2.0 = -1.1699667343192828 = 0.13275615748360492 = -0.0470241166552393 = 2.0316863985002866 = 2.0 = -1.0 State: {'0110': (0.7071067811865114+0j), '1001': (-0.7071067811865837+0j)} = 2.0 = -1.1699667343192828 = 0.13275615748360492 = -0.047024116655239305 = 2.0316863985002866 = 2.0 = 0.0 State: {'0101': (-0.3041607957254038+0j), '0110': (-0.6357045514282725+0j), '1001': (-0.6357045514282771+0j), '1010': (-0.3150327870721426+0j)} = 2.0000000000000018 = -1.1029846441120226 = 0.15484054539998238 = -0.18185168359239984 = 1.1233821369006018 = 0.0 = 0.0 State: {'0111': (0.9972637706196074+0j), '1011': (0.07392544764533501+0j)} = 2.9999999999999996 = -0.9731508943779752 = 0.1281703308946106 = -0.06264988334561404 = 2.4206463897437396 = 0.7499999999999999 = 0.5 State: {'1101': (-0.9972637706196075+0j), '1110': (-0.07392544764533537+0j)} = 3.000000000000001 = -0.9731508943779754 = 0.12817033089461063 = -0.06264988334561411 = 2.4206463897437405 = 0.7500000000000002 = -0.5000000000000001 State: {'0100': (0.960237434911092+0j), '1000': (0.2791846496413546+0j)} = 0.9999999999999997 = -0.9211167037408632 = 0.011311318681378196 = -0.05184639685022145 = 0.2422026353753374 = 0.7499999999999997 = -0.49999999999999983 State: {'0001': (-0.9602374349110921+0j), '0010': (-0.27918464964135487+0j)} = 1.0 = -0.9211167037408637 = 0.011311318681378224 = -0.051846396850221514 = 0.24220263537533723 = 0.75 = 0.5 State: {'0111': (-0.07392544764533505+0j), '1011': (0.9972637706196074+0j)} = 2.9999999999999996 = -0.8748597044244959 = 0.27009814155620415 = -0.07842246662010388 = 3.6744128057571195 = 0.7499999999999999 = 0.5 State: {'1101': (0.07392544764533537+0j), '1110': (-0.9972637706196075+0j)} = 3.000000000000001 = -0.8748597044244959 = 0.27009814155620415 = -0.07842246662010385 = 3.674412805757121 = 0.7500000000000002 = -0.5000000000000001 State: {'0101': (-0.12987265201103357+0j), '0110': (0.2593294436238296+0j '1010': (-0.9212109278548157+0j)} = 1.9999999999999996 = -0.7805602228690042 = 0.24082752532511223 = 0.011692055104520219 = 3.575907422588924 = 0.0 = 0.0 State: {'0001': (0.27918464964135487+0j), '0010': (-0.9 = 1.0 = -0.651025366468690 = 0.121 = 0.004 = 1.789 = 0.7 = 0.5 State: {'0100': (-0.27918464964135464+0 = 0.99 = -0 = 0.121 = 0.004 = 1.789 = 0.7 = -0.49 State: {'1111': (1+0j)} = 4.0 = -0 = 0.265 = -0.09 = 4.063 = 0.0 = 0.0 State: {'0000': (1+0j)} = 0.0 = 0. = 0.0 = 0.0 = 0.0 = 0.0 ```

I will focus on three states:

The ground state:

State: {'0101': (0.943726286909139+0j), '0110': (-0.169197999288364+0j), '1001': (-0.16919799928804474+0j), '
1010': (-0.22830849627297578+0j)}
        <psi|ParticleNumber|psi> = 2.000000000000001
        <psi|ElectronicEnergy|psi> = -1.1709917394922589
        <psi|DipoleMomentX|psi> = 0.0026004017257202448
        <psi|DipoleMomentY|psi> = 0.02908727852216155
        <psi|DipoleMomentZ|psi> = 1.3957696360113352
        <psi|AngularMomentum|psi> = 0.0
        <psi|Magnetization|psi> = 0.0

This has angular momentum 0 as expected. Note, its in-phase combination of the states 0110 and 1001.

The first (triple-degenerate) excited state:

State: {'0011': (1+0j)}
        <psi|ParticleNumber|psi> = 2.0
        <psi|ElectronicEnergy|psi> = -1.169966734319283
        <psi|DipoleMomentX|psi> = 0.13275615748360492
        <psi|DipoleMomentY|psi> = -0.04702411665523931
        <psi|DipoleMomentZ|psi> = 2.0316863985002866
        <psi|AngularMomentum|psi> = 2.0
        <psi|Magnetization|psi> = 1.0

State: {'1100': (1+0j)}
        <psi|ParticleNumber|psi> = 2.0
        <psi|ElectronicEnergy|psi> = -1.1699667343192828
        <psi|DipoleMomentX|psi> = 0.13275615748360492
        <psi|DipoleMomentY|psi> = -0.0470241166552393
        <psi|DipoleMomentZ|psi> = 2.0316863985002866
        <psi|AngularMomentum|psi> = 2.0
        <psi|Magnetization|psi> = -1.0

State: {'0110': (0.7071067811865114+0j), '1001': (-0.7071067811865837+0j)}
        <psi|ParticleNumber|psi> = 2.0
        <psi|ElectronicEnergy|psi> = -1.1699667343192828
        <psi|DipoleMomentX|psi> = 0.13275615748360492
        <psi|DipoleMomentY|psi> = -0.047024116655239305
        <psi|DipoleMomentZ|psi> = 2.0316863985002866
        <psi|AngularMomentum|psi> = 2.0
        <psi|Magnetization|psi> = 0.0

All of these have angular momentum 2 and the three magnetization terms 1, -1, 0. The last one, is your $\psi_1$, since it has the out-of-phase combination of the terms 0110 and 1001. Given that it is part of this triple-degenerate excited state, I hope this convinces you that the correct angular momentum for $\pi_1$ is in fact 2.

The next excited state:

State: {'0101': (-0.3041607957254038+0j), '0110': (-0.6357045514282725+0j), '1001': (-0.6357045514282771+0j),
 '1010': (-0.3150327870721426+0j)}
        <psi|ParticleNumber|psi> = 2.0000000000000018
        <psi|ElectronicEnergy|psi> = -1.1029846441120226
        <psi|DipoleMomentX|psi> = 0.15484054539998238
        <psi|DipoleMomentY|psi> = -0.18185168359239984
        <psi|DipoleMomentZ|psi> = 1.1233821369006018
        <psi|AngularMomentum|psi> = 0.0
        <psi|Magnetization|psi> = 0.0

Finally, we have another in-phase combination of 0110 and 1001 in this next excited state, which has again angular momentum 0. This comes closest to your $\psi_2$.

I hope this clears up any confusion you may have had.

MartinBeseda commented 2 years ago

@mrossinek Oh, thank you very much! It is much more clear, as I got puzzled by the plus sign instead of minus one... This one is caused by the qubit ordering adopted by Qiskit, is it right?

mrossinek commented 2 years ago

In the Jordan-Wigner mapping case, the 4 qubits correspond to:

0101
dcba

where:

a: alpha-spin orbital of MO index 0
b: alpha-spin orbital of MO index 1
c: beta-spin orbital of MO index 0
d: beta-spin orbital of MO index 1

The reason being:

I cannot say with certainty whether this affects the sign compared to a big-endian qubit register but I do find it likely due to the systems parity.

I am closing this issue as answered :+1: