search

Exa-sCI / QuantumEnvelope

7 stars 11 forks source link

PT2 reference #14

Closed kgasperich closed 4 years ago

kgasperich commented 4 years ago

@scemama is there a way to compute the PT2 deterministically in QP2? In our tests, are we comparing to the stochastic PT2 value?

scemama commented 4 years ago

The PT2 is a hybrid between deterministic and stochastic. So you can compute the exact PT2 if you let it run long enough. 2 possibilities:

  1. If you set the pt2_relative_error to zero, it will stop when the error bar is zero -> when the deterministic value is obtained.
  2. you can set do_pt2 = False to use the old deterministic algorithm. Then, you need to set threshold_generators = 1.0.

Note : We use the 2x2 diagonalization formula. You can change QP to use the standard formula in cipsi/selection.irp.f (here : https://github.com/QuantumPackage/qp2/blob/08f3ba0d54d8d72000cb973706f885d4f61f1e2b/src/cipsi/selection.irp.f#L745). This is what I did to get the reference value in the notebook. I also computed the PT2 with my Quantum Camel and the value agrees with QP. But I can't find the same PT2 energy with the miniapp.

kgasperich commented 4 years ago

I tried both of those methods in QP2, and I get two different values for the PT2. Also, neither of them agrees with either of the PT2 energies calculated with the miniapp.

maybe I made a mistake in the PT2 formula in QP2?

val = alpha_h_psi * alpha_h_psi
e_pert = val/delta_E
scemama commented 4 years ago

Here is the code I use:

      do istate=1,N_states                                                                    
        delta_E = E0(istate) - Hii + E_shift
        alpha_h_psi = mat(istate, p1, p2)
        val = alpha_h_psi + alpha_h_psi
        tmp = dsqrt(delta_E * delta_E + val * val)
        if (delta_E < 0.d0) then
            tmp = -tmp
        endif
        e_pert = 0.5d0 * (tmp - delta_E)
        if (dabs(alpha_h_psi) > 1.d-4) then
          coef = e_pert / alpha_h_psi
        else
          coef = alpha_h_psi / delta_E
        endif
    val = alpha_h_psi * alpha_h_psi   ! <----
    e_pert = val/delta_E              ! <----
        pt2(istate) = pt2(istate) + e_pert
        variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
        norm(istate) = norm(istate) + coef * coef

This is what I get for 20 determinants of F2 6-31G:

Stochastic PT2 with pt2_relative_error = 0.

 N_det             =           20
 N_states          =            1
 N_sop             =            7

 * State            1
 < S^2 >         =  -3.377316310239694E-021
 E               =   -198.709587665653     
 Variance        =   0.965569669809728     
 PT norm         =   0.284185791926475     
 PT2             =  -0.252991965069428     
 rPT2            =  -0.234086752729398     
 E+PT2            =   -198.962579630723       +/-   0.000000000000000E+000
 E+rPT2           =   -198.943674418383       +/-   0.000000000000000E+000

 -----
.. >>>>> [ IO READ: io_mo_integrals_pseudo ] <<<<< ..

.. >>>>> [ RES  MEM :       0.042171 GB ] [ VIRT MEM :       2.949604 GB ] <<<<< ..
.. >>>>> [ WALL TIME:       0.138875  s ] [ CPU  TIME:       0.980000  s ] <<<<< ..

 Energy components
 =================

 Vnn  : Nucleus-Nucleus   potential energy
 Ven  : Electron-Nucleus  potential energy
 Vee  : Electron-Electron potential energy
 Vecp : Potential energy of the pseudo-potentials
 T    : Electronic kinetic energy

 State            1
 ---------

 Vnn  =    30.3586331075289     
 Ven  =   -537.257053199709     
 Vee  =    109.388886965744     
 Vecp =   0.000000000000000E+000
 T    =    198.799945460782

Deterministic PT2 with threshold_generators = 1.0:

 N_det             =           20
 N_states          =            1
 N_sop             =            7

 * State            1
 < S^2 >         =  -2.017835100698860E-023
 E               =   -198.709587665565     
 Variance        =   0.965569669809790     
 PT norm         =   0.284185791937197     
 PT2             =  -0.252991965076561     
 rPT2            =  -0.234086752734678     
 E+PT2 (approx)   =   -198.962579630642       +/-   1.422909060022790E-321
 E+rPT2(approx)   =   -198.943674418300       +/-   0.000000000000000E+000

 -----
.. >>>>> [ IO READ: io_mo_integrals_pseudo ] <<<<< ..

.. >>>>> [ RES  MEM :       0.041672 GB ] [ VIRT MEM :       2.952419 GB ] <<<<< ..
.. >>>>> [ WALL TIME:       0.139173  s ] [ CPU  TIME:       0.972000  s ] <<<<< ..

 Energy components
 =================

 Vnn  : Nucleus-Nucleus   potential energy
 Ven  : Electron-Nucleus  potential energy
 Vee  : Electron-Electron potential energy
 Vecp : Potential energy of the pseudo-potentials
 T    : Electronic kinetic energy

 State            1
 ---------

 Vnn  =    30.3586331075289     
 Ven  =   -537.257053199503     
 Vee  =    109.388886965703     
 Vecp =   0.000000000000000E+000
 T    =    198.799945460706

The differences come from round-off errors.

Here is the input if you want to try (SCF MOs):

  2

  F             0.00000000       0.00000000       0.00000000
  F             0.00000000       0.00000000       1.41190000

E(SCF) = 198.6460967431443

  0.967460412647
    +++++++++-------------------------------------------------------
    +++++++++-------------------------------------------------------

  -0.251917601315
    ++++++-+++------------------------------------------------------
    ++++++-+++------------------------------------------------------

  0.0149434506023
    +++++++++-------------------------------------------------------
    ++++++-++-------+-----------------------------------------------

  0.0149434506023
    ++++++-++-------+-----------------------------------------------
    +++++++++-------------------------------------------------------

  0.00389378158862
    +++++++-+--+----------------------------------------------------
    ++++++-+++------------------------------------------------------

  0.00389378158862
    ++++++-+++------------------------------------------------------
    +++++++-+--+----------------------------------------------------

  -0.0038692589098
    ++++++-+++------------------------------------------------------
    ++++++++--+-----------------------------------------------------

  -0.0038692589098
    ++++++++--+-----------------------------------------------------
    ++++++-+++------------------------------------------------------

  -0.00354439080677
    ++++++--++-+----------------------------------------------------
    +++++++++-------------------------------------------------------

  -0.00354439080677
    +++++++++-------------------------------------------------------
    ++++++--++-+----------------------------------------------------

  0.00352112632768
    ++++++-+-++-----------------------------------------------------
    +++++++++-------------------------------------------------------

  0.00352112632768
    +++++++++-------------------------------------------------------
    ++++++-+-++-----------------------------------------------------

  0.000416235761727
    ++++++-++----+--------------------------------------------------
    +++++++++-------------------------------------------------------

  0.000416235761722
    +++++++++-------------------------------------------------------
    ++++++-++----+--------------------------------------------------

  0.000349390781845
    +++++++-++------------------------------------------------------
    ++++++-++--+----------------------------------------------------

  0.000349390781843
    ++++++-++--+----------------------------------------------------
    +++++++-++------------------------------------------------------

  -0.000348132582123
    ++++++++-+------------------------------------------------------
    ++++++-++-+-----------------------------------------------------

  -0.000348132582118
    ++++++-++-+-----------------------------------------------------
    ++++++++-+------------------------------------------------------

  0.000345078832927
    +++++++-+------+------------------------------------------------
    +++++++++-------------------------------------------------------

  0.000345078832914
    +++++++++-------------------------------------------------------
    +++++++-+------+------------------------------------------------
TApplencourt commented 4 years ago

I also computed the PT2 with my Quantum Camel and the value agrees with QP. But I can't find the same PT2 energy with the miniapp.

Oh, that means we have an error in the miniapp for sure. This a progress

scemama commented 4 years ago

@TApplencourt Yes, it is a bug in the miniapp. It is probably in the generation of all singles and doubles, avoiding double counting and excluding those that are already in the internal space. If the variational energy is OK on a large CI, it means that the Slater rules are correct.

TApplencourt commented 4 years ago

Let me put more unit test <3

TApplencourt commented 4 years ago

Solved. It was an iterator issue. Everything is working know

scemama commented 4 years ago

PT2 is OK?

On 03/03/2020 20:29, Thomas Applencourt wrote:

Solved. It was an iterator issue. Everything is working know

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/TApplencourt/QuantumEnvelope/issues/14?email_source=notifications&email_token=ABNRVYULADKNFHNME2XR2JLRFVLDDA5CNFSM4K76GXV2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOENU2MXY#issuecomment-594126431, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABNRVYV3JGDRMSC3XBMVWXTRFVLDDANCNFSM4K76GXVQ.

TApplencourt commented 4 years ago

For one example (the 10 determinant) it's working. But it doesn't seem to work for others. the reference PT2 energy is wrong? Can you do more check on your end?

On Tue, Mar 3, 2020 at 2:58 PM Anthony Scemama notifications@github.com wrote:

PT2 is OK?

On 03/03/2020 20:29, Thomas Applencourt wrote:

Solved. It was an iterator issue. Everything is working know

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub < https://github.com/TApplencourt/QuantumEnvelope/issues/14?email_source=notifications&email_token=ABNRVYULADKNFHNME2XR2JLRFVLDDA5CNFSM4K76GXV2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOENU2MXY#issuecomment-594126431 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/ABNRVYV3JGDRMSC3XBMVWXTRFVLDDANCNFSM4K76GXVQ .

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/TApplencourt/QuantumEnvelope/issues/14?email_source=notifications&email_token=ABRY72ZZXXSLEFWSDISQAQ3RFVVPXA5CNFSM4K76GXV2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOENVD3FA#issuecomment-594165140, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABRY72YAEM6OAXDFYP2OABTRFVVPXANCNFSM4K76GXVQ .

scemama commented 4 years ago

Done !