Multi-geminal wave functions

[Vladimir_Konjkov]

Patch now available in current beta.

Best,
Pablo

Hello Pablo, thanks for the opportunity to test the calculations with geminals ansatz.

I chose the simplest systems and performed the calculations in the qchem program.

1. Li-atom in cc-pVDZ basis, WFN consists of 1 geminal orbitals and an orbital filled with an unpaired electron, so QCHEM input/output and CASINO input/output in attachment.

Final energy = -7.432605555957891

Energy decomposition, spin= TRUE

Nuclear Attr. Energy = 0.0000000000

One-El. Pot. Energy = -17.1450277592
One-El. Kin. Energy = 7.4317293525
Tot. Intra-gem Energy = -9.2971415985
Inter-gem Coul. Energy = 2.6905114013
Inter-gem Exch. Energy = -0.8259753588
Correlation Energy = -0.0004001996
VAR verifying total -7.432605555957891 error: -0.000000000000001

Geminal occupations, energies and coefficients.

Geminal 1 E = -3.387669043962796
2.00[s],
0.99996729 -0.00464321 -0.00464307 -0.00464305 -0.00048555 -0.00031930 -0.00031926 -0.00031921 -0.00031917 -0.00031659 -0.00001423 -0.00001403 0.00001401.

Occupations and energies of unpaired orbitals.
1.00[s],
-2.18040047

Unfortunately CASINO cannot accept an unpaired electron on second orbital and returns an error: Kinetic energy test failed: analytical derivatives misbehave.

2. Be-atom in cc-pVDZ basis, WFN consists of 2 geminals strongly ortogonal, so QCHEM input/output and CASINO input/output in attachment.

Final energy = -14.617016246569332

Energy decomposition, spin= TRUE

Nuclear Attr. Energy = 0.0000000000

One-El. Pot. Energy = -33.6940528313
One-El. Kin. Energy = 14.6150473252
Tot. Intra-gem Energy = -16.5188194534
Inter-gem Coul. Energy = 1.9430696836
Inter-gem Exch. Energy = -0.0412664768
Correlation Energy = -0.1214694892
VAR verifying total -14.617016246569332 error: 0.000000000000017

Geminal occupations, energies and coefficients.

Geminal 1 E = -11.707720511899904
2.00[s],
0.99999849 -0.00100242 -0.00100221 -0.00100180.

Geminal 2 E = -1.007492527867571
1.82[s],
0.95271649 -0.17380094 -0.17379549 -0.17379340 -0.03715659 -0.00818632 -0.00818536 -0.00818511 -0.00818199 -0.00818198

each virtual orbital belongs either to the first or second geminal, this results in orthogonality, and I manually distributed virtual orbitals based on the QCHEM output file, however, the energy does not coincide with that found in QCHEM output file.
I will be very grateful if you describe the algorithm for constructing a wave function (GEMINAL section) for two orthogonal geminal.

if changes were made to a separate git branch it would be much easier to test, without breaking master branch.

Best Vladimir

[Pablo_Lopez_Rios]

Hi Vladimir,

Li-atom in cc-pVDZ basis, WFN consists of 1 geminal orbitals and an orbital filled with an unpaired electron

Unpaired electrons are indeed not supported; the geminals are supposed to be determinants of square matrices of size Nup x Ndown which must therefore satisfy Nup = Ndown. If you could point at a description of how geminals are defined for systems with Nup /= Ndown I might be able to implement this.

Be-atom in cc-pVDZ basis, WFN consists of 2 geminals strongly ortogonal,

You could start by trying each of these geminals separately to see if individual energies match those from QCHEM. In the CASINO input files you have set the coefficients of both geminals to one, but I doubt this is optimal. Does QCHEM provide the coefficients? If not, it should be easy to fix c_1=1 and optimize c_2 within VMC (having first verified the energies of the individual geminals).

if changes were made to a separate git branch it would be much easier to test, without breaking master branch.

CASINO only uses a master branch; this is out of my control. In any case, nothing has been broken in the master branch (my commits pass the autotest suite).

Best,
Pablo

[Pablo_Lopez_Rios]

Hi Vladimir,

Forgot to reply to this:

I will be very grateful if you describe the algorithm for constructing a wave function (GEMINAL section) for two orthogonal geminal.

I have no experience with constructing multigeminal wave functions. I would assume that https://www.tcm.phy.cam.ac.uk/~pl275/do ... thesis.pdf might shed some light on this (for the specific case of the electron gas)?

Best,
Pablo

[Vladimir_Konjkov]

Hi Vladimir,

Li-atom in cc-pVDZ basis, WFN consists of 1 geminal orbitals and an orbital filled with an unpaired electron

Unpaired electrons are indeed not supported; the geminals are supposed to be determinants of square matrices of size Nup x Ndown which must therefore satisfy Nup = Ndown. If you could point at a description of how geminals are defined for systems with Nup /= Ndown I might be able to implement this.

QCHEM uses antisymmetrized product of strongly orthogonal geminals (APSG) ansatz introdused by Rassolov in https://aip.scitation.org/doi/10.1063/1.1503773

To summarize, the complete specification of the APSSG (or SSG) model for the system of n_alpha electrons with spin up and n_beta electrons with the spin down (we assume n_alpha > n_beta ) is as follows:
The wave function has the form (5):

psi SSG = A[geminal_1(r1, r2) ... geminal_n_beta(r2n_beta-1, r2n_beta) * phi_i(r2n_beta+1)...phi_j(rn_beta + n_alpha)]

where A is an operator that antisymmetrizes the product in square brackets with respect to all electron permutations i.e. determinant.
The difference from your case is that singly occupied orbitals included as ordinary single-electron spin-orbitals - phi().

Geminals in turn can be represented as:

geminal(r1, r2) = sum Dij * A[phi(r1)phi(r2)]

where phi is spin-orbitals and i < j restriction is introduced in order to prevent the double-counting of configurations.
This equation can be re-written with the spin and spatial functions separated for a singlet geminal in the MO basis:

geminal(r1, r2) = sum Dij * F_i(r1)*F_j(r2)*(alpha(s1)*beta(s2)-beta(s2)*alpha(s1))

Spin function is antisymmetric, and so the coefficients Dij must be symmetric in order to preserve the whole wavefunction’s antisymmetry. This allows the spatial function to be reduced to diagonal form.

geminal(r1, r2) = sum Dk * PNO_k(r1)*PNO_k(r2)*(alpha(s1)*beta(s2)-beta(s2)*alpha(s1))

where PNO is (pseudo)-natural orbitals and k belongs to Arai subspaces - unintersecting supspaces of MOs which enforce strong orthogonality of the geminals.

Imposing the strong orthogonality constraint means that the one-electron space spanned by PNO is factored into n disjoint subspaces, such that each geminal is expanded in the one-electron functions belonging to that geminal’s subspace.
These subspaces are termed Arai’s subspaces, after Arai’s theorem which outlines this concept. In order for the wavefunction to be adapted to open-shell systems, the geminals are multiplied by one electron orbitals;
The one-electron orbitals are optimized under a further orthogonalization constraint;

As far as I understand QCHEM provides PNOs orbitals in MOLDEN file and Dk for every geminals as a list in output, for Be example 1-st geminal consists of 1, 12, 13, 14 orbital and 2-nd geminal consist of from 2 to 11 orbitals, I copied D weights into GEMINAL section of casl file from QCHEM output, that is, in the calculation without JASTROW factor I have to get the same energy as QCHEM output. I could not find weight of each geminal but I think they are equal to one.

I would assume that https://www.tcm.phy.cam.ac.uk/~pl275/do ... thesis.pdf might shed some light on this (for the specific case of the electron gas)?

thanks, I started to read, it is interesting, I didn’t know that there is a connection between geminals and the BCS theory, can they be used to describe Mott insulators?

Best Vladimir.

[Vladimir_Konjkov]

Hello Pablo,

Returning to this thread after a few years. I have been testing the MAGP wave function again (CASINO v2.14.0, current beta) on the Be atom with Gaussian orbitals, and I believe I found a remaining bug in geminal.f90: stale one-electron log-gradients (FIGEM) served from scratch buffers after accepted single-electron moves.

Symptom

Be atom, all-electron, cc-pVQZ, a simple diagonal GEMINAL block, plain VMC. The three kinetic-energy estimators disagree far beyond statistics (they must agree for exact wave-function derivatives):

Code: Select all

without Jastrow:  KEI = 14.61832(45)   TI = 14.62580(31)   FISQ = 14.63328(52)
with Jastrow:     KEI = 14.7899        TI = 14.7955        FISQ = 14.8011

FISQ - KEI is ~15 mHa at ~20 sigma in both cases, with both vmc_method 1 and 3. The built-in kinetic-energy check at VMC startup does not catch it because the problem only appears via the accepted-move buffer-copy path.

Cause

Three related indexing problems:

1) accept_move_geminal, block "Update FIGEM and LAPGEM": the validity flags are copied for ALL electrons, but the data are copied only for the moved electron - and addressed with the spin-relative index ie, while figem_scr/figem_valid are dimensioned over the absolute electron index (3,ngems,netot,nscratch). After an accepted move the destination buffer claims valid per-electron gradients that it does not actually hold (or holds in the wrong slot for down-spin electrons).

2), 3) get_figem and get_lapgem: the guard before refreshing the geminal inverse tests figem_valid(igem,ie,is) with the spin-relative ie where the
absolute ii is meant, so for down-spin electrons the wrong electron's flag is consulted and the inverse may not be recomputed when it must be.

Fix

Code: Select all

--- a/src/geminal.f90
+++ b/src/geminal.f90
@@ -1249,9 +1249,11 @@

 ! Update FIGEM and LAPGEM.
  figem_valid(:,:,is)=figem_valid(:,:,js)
- do igem=1,ngems
-  if(figem_valid(igem,ie,js))figem_scr(:,igem,ie,is)=figem_scr(:,igem,ie,js)
- enddo ! igem
+ do i=1,netot
+  do igem=1,ngems
+   if(figem_valid(igem,i,js))figem_scr(:,igem,i,is)=figem_scr(:,igem,i,js)
+  enddo ! igem
+ enddo ! i

  if(.not.use_backflow)then
   lapgem_valid(:,:,is)=lapgem_valid(:,:,js)
@@ -2386,7 +2388,7 @@
  ie=which_ie(ii) ; ispin=which_spin(ii)

  do igem=1,ngems
-  if(.not.figem_valid(igem,ie,is))then
+  if(.not.figem_valid(igem,ii,is))then
    call get_gem_inv(igem,is)
   endif
  enddo
@@ -2567,7 +2569,7 @@
  ie=which_ie(ii) ; ispin=which_spin(ii)

  do igem=1,ngems
-  if(.not.figem_valid(igem,ie,is))then
+  if(.not.figem_valid(igem,ii,is))then
    call get_gem_inv(igem,is)
   endif
  enddo

Verification

- After the patch the KEI/TI/FISQ triplet agrees within error bars on the same tests (with and without Jastrow, vmc_method 1 and 3).

- The wave-function values themselves were never wrong: a line scan (runtype : plot) matches an independent numpy implementation of
sum_n c_n det[Phi_n] up to a constant normalization at every point.

- An end-to-end physics check: for Be the 3-geminal expansion

Code: Select all

G1 = phi_1s x phi_1s + phi_2s x phi_2s            (c = +1)
G2 = phi_1s x phi_1s + eps * sum_p phi_2p x phi_2p (c = +1)
G3 =                   eps * sum_p phi_2p x phi_2p (c = -1)

is algebraically identical to the standard 4-determinant 2s^2 -> 2p^2 CSF expansion (the third geminal cancels the spurious eps^2 double-promotion
pair products of the second, which otherwise cost ~3*eps^4*dE with an empty 1s core). With the patch, VMC gives E = -14.617(2) for the geminal form vs
E = -14.6168(6) for the equivalent MDET run with the same CASSCF natural orbitals - exact agreement.

Best regards,
Vladimir.

[Vladimir_Konjkov]

A follow-up on the constrained optimization of this 3-geminal ansatz, plus one more (cosmetic) issue found on the way.

Optimizing with the g elements tied together works nicely:

Code: Select all

Constraints:
  Equate 1: [ Diagonal: 3:5, Geminals: 2:3 ]

with exactly one element of the constrained set flagged "optimizable" and the remaining five flagged "determined" (as check_g_constraint requires). The
attached archive contains the complete run directory (input, gwfn.data, parameters.casl with its iterations, out) for Be with a Jastrow factor and emin.

The cosmetic issue: in the parameters.N.casl files written out during optimization, only the reference element of the constraint is updated. E.g. after a few cycles I get

Code: Select all

Geminal 2:
  Parameters:
    c: [ 1.0000000000000000, fixed ]
    g_1,1: [ 1.0000000000000000, fixed ]
    g_3,3: [ -0.17922936592613789, optimizable ]
    g_4,4: [ -0.18970000000000001, determined ]
    g_5,5: [ -0.18970000000000001, determined ]

where the "determined" entries keep their initial value -0.1897 instead of following the reference. This is only in the written file: in memory apply_constraints is called after every parameter update, and read_geminal also calls apply_constraints right after parsing, so both the running wave function and any restart from such a file are correct - the file is just misleading to read.

The reason is in update_geminal_casl, which skips both fixed and determined elements when updating the g blocks:

Code: Select all

if(gmat_opt(irow,icol,igem)==opt_fixed.or.&
 &gmat_opt(irow,icol,igem)==opt_determined)cycle

while the c-coefficient block right above it DOES write determined values (its condition includes opt_determined). Making the two consistent - e.g. cycling only on opt_fixed in the g loop - would make the written casl reflect the actual parameter values.

Best regards,
Vladimir.

[Vladimir_Konjkov]

One more follow-up, closing the Be story: the MAGP wave function can build the correct multi-determinant nodal surface entirely on its own, starting from HF orbitals - no CASSCF natural orbitals in gwfn.data needed.

The idea: for a symmetric pairing matrix, g = U lambda U^T, so the off-diagonal elements of g ARE orbital rotations. Freeing an off-diagonal block should therefore do the job of CASSCF orbital optimization inside VMC.

Setup

HF/cc-pVQZ gwfn.data for Be. The p-type MOs come in 4 radial shells per component: x = (5, 9, 17, 42), y = (4, 8, 15, 40), z = (3, 7, 16, 41).
Geminal 2 carries 1s x 1s plus a full symmetric 4x4 block over the p MOs of each component; Geminal 3 (c = -1) mirrors the p blocks to cancel the spurious both-pairs-in-p configurations discussed above. The 10 unique parameters (upper triangle of the x block) are the only free ones; long-form constraints tie each to its transpose, to the y/z copies (S symmetry) and to the Geminal-3 mirror, e.g.

Code: Select all

Constraints:
  2^g_5,5=2^g_4,4=2^g_3,3=3^g_5,5=3^g_4,4=3^g_3,3
  2^g_5,9=2^g_9,5=2^g_4,8=2^g_8,4=2^g_3,7=2^g_7,3=3^g_5,9=3^g_9,5=3^g_4,8=3^g_8,4=3^g_3,7=3^g_7,3
  ...

The full run directory is attached.

Results (all-electron VMC, Jastrow, 10^6 steps)

varmin systematically stalls at E ~ -14.65: the variance is almost insensitive to the near-degeneracy weight and to the radial shape of the correlating 2p, a well-known failure mode for CSF-type weights. Energy minimization does the job:

Code: Select all

VMC #1: E = -14.575(2)    var = 2.9(2)     start (Jastrow off)
VMC #2: E = -14.6514(3)   var = 0.0370(3)  varmin, Jastrow only
VMC #3: E = -14.6663(2)   var = 0.0291(4)  emin cycle 1
VMC #4: E = -14.6666(2)   var = 0.0281(2)  emin cycle 2
4-det MDET control (CASSCF orbitals + Jastrow): E = -14.6669(9)

Statistically identical to the multi-determinant control.

What the optimizer built

Eigendecomposition of the converged 4x4 p block:

Code: Select all

lambda_1 = -0.1682   eigenvector (2p,3p,4p,5p)_HF = (-0.83, -0.55, +0.11, 0.00)
lambda_2 = -0.0066   (|l2/l1| = 3.9%)
lambda_3,4 ~ 3e-4

The block is essentially rank 1, and the dominant eigenvector is the compact correlating 2p: 83% HF-2p reshaped with a 55% admixture of HF-3p -
exactly the radial contraction that CASSCF orbital optimization would perform, here done by the off-diagonal g elements during VMC optimization.
The tiny second channel is a bonus beyond CAS(2,4).

So in practice: with the constraint machinery, MAGP + emin recovers the multi-determinant nodal surface from an HF starting point, with the orbital
optimization built into the g matrix. varmin should be avoided for the geminal parameters themselves.

Best regards,
Vladimir.

[Vladimir_Konjkov]

A second, unrelated bug: MAGP crashes with an FPE for any decently-sized atom, even with a plain diagonal (HF-equivalent) GEMINAL block. Found on Ne (N=5, all-electron cc-pVQZ); reproduces on a clean, unpatched geminal.f90 (orthogonal to the FIGEM fix above).

Cause

Code: Select all

REAL(dp),PARAMETER :: tol_log_softzero=-30._dp

used in caldet to flag a geminal determinant as "hard zero" (singular). The equivalent constant in slater.f90 is -690._dp - 23 orders of magnitude stricter in log-space. For a 5-electron atom with a reasonably large basis (cc-pVQZ), the perfectly valid, non-singular 5x5 geminal determinant routinely computes to log(det) ~ -33 to -37 (nowhere near underflow), tripping the -30 cutoff but not Slater's -690 one. The code then discards the only geminal at nearly every sampled configuration, collapsing the wave function to zero - symptom: ke_verbose reports "Analytical derivatives misbehave", then FPE abort on VMC start.

Fix

Code: Select all

- REAL(dp),PARAMETER :: tol_log_softzero=-30._dp
+ REAL(dp),PARAMETER :: tol_log_softzero=-690._dp

Verification

ke_verbose now reports "Geminals - gradient: optimal, Laplacian: optimal". VMC energy for the diagonal (HF-equivalent) geminal, no Jastrow: E = -128.545(2), matching the ORCA HF/cc-pVQZ reference -128.5434696 within 1 sigma.

Best regards,
Vladimir.

[Vladimir_Konjkov]

A third bug, MAGP + backflow only: heap corruption from a wrongly dimensioned validity array.

Symptom

Ne, all-electron cc-pVQZ, psi_s : geminal, backflow : T, varmin. The optimization completes and writes correlation.out, but the job aborts on exit with "free(): invalid next size (fast)" and a backtrace ending in __libc_free from main - glibc only notices the damage when the arrays are deallocated. A make debug binary points straight at the source:

Code: Select all

  At line 4397 of file geminal.f90 
  Fortran runtime error: Index '2' of dimension 2 of array 'orb_sderiv_valid'  above upper bound of 1
  

Cause

In setup_geminal, inside the if(use_backflow) block, ORB_SDERIV_VALID gets real1_complex2 as its second dimension, but it is indexed by spin everywhere it is used - in get_orbvals (line 4397, orb_sderiv_valid(i,spin,is)) and in the "Update orb_sderiv" block of accept_move_geminal (do j=1,2). Its siblings are dimensioned accordingly - orb_val_valid(nemax,2,nscratch), same for orb_grad_valid and orb_lap_valid - and ORB_SDERIV_SCR itself correctly carries nspin. For a real wave function real1_complex2 = 1, so the array is half the size it needs and every down-spin flag is written past the end of it.

Fix

Code: Select all

  --- a/src/geminal.f90                                                
  +++ b/src/geminal.f90                                                
  @@ -710,7 +710,7 @@
  
   ! orb_sderiv
     allocate(orb_sderiv_scr(6,norb,nemax,real1_complex2,nspin,nscratch),&
  -   &orb_sderiv_valid(nemax,real1_complex2,nscratch),stat=ialloc)    
  +   &orb_sderiv_valid(nemax,nspin,nscratch),stat=ialloc)
     call check_alloc(ialloc,"SETUP_GEMINAL","orb_sderiv")
  

The other flag arrays in the same block (gem_sderiv_pair_valid, Hgem_valid, Farray_valid, Harray_valid) are fine.

Verification

With the patch the bounds-checked binary runs the same Ne backflow varmin to the end with no out-of-bounds report, and the opt binary no longer aborts in free().

Best regards,
Vladimir.

[Vladimir_Konjkov]

A further bug in geminal.f90, again MAGP + backflow only (psi_s : geminal with backflow : T), on top of the ORB_SDERIV_VALID allocation fix above. Found on Ne, all-electron cc-pVQZ, varmin. Still verifying, so posting preliminarily.

Symptom

VMC converges to energies below the exact non-relativistic ground state:
E = -128.53(8) -> -128.98(2) -> -129.09(3) for Ne, against -128.9376. The kinetic-energy check reports

Code: Select all

  Geminals - gradient: optimal, Laplacian: poor.

  KEI  = 112.51(92)
  TI   = 122.34(1.49)
  FISQ = 132.17(2.85)
  

KEI and FISQ disagree by ~20 au, and KEI is what enters the total energy - so the "energy" is just a broken estimator. The Laplacian is graded "optimal" while the backflow parameters are still zero, and degrades to "good", then "poor", as they grow during the optimization.

Cause

get_Hgems_down_ii addresses the pair second derivatives with idir where kdir is meant. kdir is the running index mapping a direction pair onto the sderiv component ordering (xx, yy, zz, xy, xz, yz); the up-spin routine uses it correctly, while idir runs over {1,2,3,1,1,2}. So the xy, xz and yz components of the down-spin diagonal Hgem block read the xx, xx and yy components instead. The diagonal components coincide by accident, which is why the Laplacian is exact at zero backflow, where only the trace of Hgem is needed: once backflow is on, loglap_bf contracts Hgem with the backflow Jacobian dx/dr and the off-diagonal (alpha /= beta) components start contributing garbage - for down-spin electrons only. Spin-polarised systems (ned = 0) are unaffected.

Best regards,
Vladimir.

PS maybe

Fix

Code: Select all

  --- a/src/geminal.f90                                                
  +++ b/src/geminal.f90
  @@ -2970,14 +2970,14 @@
       idir=tdir
      endif
      tmpr=ddot(nemax,gem_inv_scr(ido,1,1,igem,is),nemax,&
  -    &gem_rsderiv_pair_scr(idir,1,ido,1,igem,is),6)                  
  +    &gem_rsderiv_pair_scr(kdir,1,ido,1,igem,is),6)
      if(complex_wf)then
       tmpr=tmpr-ddot(nemax,gem_inv_scr(ido,1,2,igem,is),nemax,&
  -     &gem_rsderiv_pair_scr(idir,1,ido,2,igem,is),6)                 
  +     &gem_rsderiv_pair_scr(kdir,1,ido,2,igem,is),6)
       tmpi=ddot(nemax,gem_inv_scr(ido,1,1,igem,is),nemax,&
  -     &gem_rsderiv_pair_scr(idir,1,ido,2,igem,is),6)+&               
  +     &gem_rsderiv_pair_scr(kdir,1,ido,2,igem,is),6)+&               
        &ddot(nemax,gem_inv_scr(ido,1,2,igem,is),nemax,&               
  -     &gem_rsderiv_pair_scr(idir,1,ido,1,igem,is),6)                 
  +     &gem_rsderiv_pair_scr(kdir,1,ido,1,igem,is),6)                 
       Hgem_scr(idir,jdir,2,ii,ii,igem,is)=tmpi                        
      endif
      Hgem_scr(idir,jdir,1,ii,ii,igem,is)=tmpr