We present an efficient approach to the electron correlation problem that is well-suited for strongly interacting many-body systems, but requires only mean-field-like computational cost. The performance of our approach is illustrated for the one-dimensional Hubbard model for different ring lengths, and for the non-relativistic quantum chemical Hamiltonian exploring the symmetric dissociation of the H50 hydrogen chain.The accurate description of the electron-electron interaction at the quantum-mechanical level is a key problem in condensed matter physics and quantum chemistry. Since most of the quantum many-body problems are extraordinarily difficult to solve exactly, different approximation schemes emerged [5][6][7][8][9], among which the density matrix renormalization group (DMRG) algorithm [10][11][12] gained a lot of popularity in both condensed matter physics [11] and quantum chemistry [13][14][15][16][17][18][19][20] over the last decade. Since the DMRG algorithm optimizes a matrix product state wavefunction, it is optimally suited for one-dimensional systems; though DMRG studies on higher-dimensional and compact systems have been reported [14,17,19,21]. Yet, novel theoretical approaches are desirable that can accurately describe strong correlation effects between electrons where the dimension of the Hilbert space exceeds the present-day limit of DMRG or general tensor-network approaches [22] allowing approximately 100 sites or 60 (spatial) orbitals, respectively.Another promising approach, suitable for larger strongly-correlated electronic systems, uses geminals (two-electron basis functions) as building blocks for the wavefunction [23][24][25][26][27][28][29][30]. One of the simplest practical geminal approaches is the antisymmetric product of 1-reference-orbital geminals (AP1roG) [31][32][33]. Unique among geminal methods, AP1roG can be rewritten as a fully general pair-coupled-cluster doubles wavefunc- * ayers@mcmaster.ca † dimitri.vanneck@ugent.be tion [34], i.e.where a † pσ and a pσ (σ = ↓, ↑) are the fermionic creation and annihilation operators, and |Φ 0 is some independent-particle wavefunction (usually the HartreeFock determinant). Indices i and a correspond to virtual and occupied sites (orbitals) with respect to |Φ 0 , P and K denote the number of electron pairs (P = N/2 with N being the total number of electrons) and orbitals, respectively, and {c a i } are the geminal coefficients. This wavefunction ansatz is size-extensive and has mean-field scaling, O P 2 (K − P ) 2 for the projected Schrödinger equation approach [31].To ensure size-consistency, however, it is necessary to optimize the one-electron basis functions [31], where all non-redundant orbital rotations span the occupiedoccupied, occupied-virtual, and virtual-virtual blocks with respect to the reference Slater determinant |Φ 0 . We have implemented a quadratically convergent algorithm: we minimize the energy with respect to the choice of the one-particle basis functions, subject to the constraint that the projected Schrödinger equati...