Automatic Routing of Goldstone Diagrams using Genetic Algorithms

TL;DR

This paper introduces a genetic algorithm-based method to automatically route complex Goldstone diagrams, optimizing their graphical representation by minimizing a cost function, and demonstrates its effectiveness on advanced quantum chemistry diagrams.

Contribution

It presents a novel genetic algorithm approach for automatic routing of Goldstone diagrams, handling both discrete and continuous parameters to produce diagrams comparable to manual routing.

Findings

Successfully routes complex diagrams with multiple loops and crossings.

Produces diagrams qualitatively equivalent to manual routing.

Applicable to various Coupled Cluster methods and perturbative expansions.

Abstract

This paper presents an algorithm for an automatic transformation (=routing) of time ordered topologies of Goldstone diagrams (i.e. Wick contractions) into graphical representations of these topologies. Since there is no hard criterion for an optimal routing, the proposed algorithm minimizes an empirically chosen cost function over a set of parameters. Some of the latter are naturally of discrete type (e.g. interchange of particle/hole lines due to antisymmetry) while others (e.g. x,y-position of nodes) are naturally continuous. In order to arrive at a manageable optimization problem the position space is artificially discretized. In terms of the (i) cost function, (ii) the discrete vertex placement, (iii) the interchange of particle/hole lines the routing problem is now well defined and fully discrete. However, it shows an exponential complexity with the number of vertices suggesting to apply a genetic algorithm for its solution. The presented algorithm is capable of routing non trivial (several loops and crossings) Goldstone diagrams. The resulting diagrams are qualitatively fully equivalent to manually routed ones. The proposed algorithm is successfully applied to several Coupled Cluster approaches and a perturbative (fixpoint iterative) CCSD expansion with repeated diagram substitution.

Figures

Tables

Table 1. Table 1: Static costs of Goldstone diagrams during gene evaluation.

Index	Example	Costs
(A)		$-10$
(B)		$+5$
(C)		$+10$
(D)		$+2\cdot 1.5$
(E)		$+100$
(F)		$-1$

Table 2. Table 2: Possible Goldstone diagram mutations during phase II.

(A)		$\leftrightarrow$
(B)		$\leftrightarrow$
(C)		$\leftrightarrow$
(D)		$\leftrightarrow$

Table 3. Table 3: Collection of optimized Goldstone diagrams for different CC approaches. For each contraction pattern only one participating diagram is shown.

Bra	Ket	Diagram
$⟨{\Psi }_{IJK}^{ABC}|$	$\text{contraction}\left[1ex\right]{\hat{V}}_N{\hat{T}}_1\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{T}}_1{\hat{T}}_2\text{contraction}\left[2ex\right]{\hat{V}}_N{\hat{T}}_1{\hat{T}}_2\text{contraction}\left[2.5ex\right]{\hat{V}}_N{\hat{T}}_1{\hat{T}}_2{\hat{T}}_2$ ${\hat{V}}_N{\hat{T}}_1{\hat{T}}_2{\hat{T}}_2{\Psi }_0⟩$111standard CC
$⟨{\Psi }_{IJKL}^{ABCD}|$	$\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{T}}_2\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{T}}_2{\hat{T}}_2$ ${\hat{V}}_N{\hat{T}}_2{\hat{T}}_2{\Psi }_0⟩$111standard CC
$⟨{\Psi }_0|$	$\text{contraction}\left[1.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N\text{contraction}\left[1.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N\text{contraction}\left[2.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N{\hat{T}}_1\text{contraction}\left[2.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N{\hat{T}}_1\text{contraction}\left[3.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N{\hat{T}}_1{\hat{T}}_2\text{contraction}\left[3.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N{\hat{T}}_1{\hat{T}}_2\text{contraction}\left[4.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N{\hat{T}}_1{\hat{T}}_2\text{contraction}\left[4.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N{\hat{T}}_1{\hat{T}}_2$ ${\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{F}}_N{\hat{T}}_1{\hat{T}}_2{\Psi }_0⟩$222expectation value/variational CC
$⟨{\Psi }_0|$	$\text{contraction}\left[1.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2\text{contraction}\left[1.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2\text{contraction}\left[2.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N\text{contraction}\left[2.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N\text{contraction}\left[3.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2\text{contraction}\left[3.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2\text{contraction}\left[4.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2{\hat{T}}_2\text{contraction}\left[4.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2{\hat{T}}_2\text{contraction}\left[5.0ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2{\hat{T}}_2\text{contraction}\left[5.5ex\right]{\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2{\hat{T}}_2$ ${\hat{T}}_2^{†}{\hat{T}}_2^{†}{\hat{V}}_N{\hat{T}}_2{\hat{T}}_2{\Psi }_0⟩$222expectation value/variational CC
$⟨{\Psi }_0|$	$\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{T}}_1\text{contraction}\left[2.0ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{T}}_1$ ${\hat{V}}_N{\hat{F}}_N{\hat{T}}_1{\Psi }_0⟩$333variance calculation
$⟨{\Psi }_0|$	$\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{V}}_N\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{V}}_N\text{contraction}\left[2.0ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{T}}_1\text{contraction}\left[2.5ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{T}}_1{\hat{T}}_1\text{contraction}\left[3.0ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{T}}_1{\hat{T}}_1{\hat{T}}_1\text{contraction}\left[3.5ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{T}}_1{\hat{T}}_1{\hat{T}}_1$ ${\hat{V}}_N{\hat{V}}_N{\hat{T}}_1{\hat{T}}_1{\hat{T}}_1{\Psi }_0⟩$333variance calculation
$⟨{\Psi }_{IJ}^{AB}|$	$\text{contraction}\left[1.0ex\right]{\hat{W}}_N{\hat{T}}_2\text{contraction}\left[1.5ex\right]{\hat{W}}_N{\hat{T}}_2\text{contraction}\left[2.0ex\right]{\hat{W}}_N{\hat{T}}_2$ ${\hat{W}}_N{\hat{T}}_2{\Psi }_0⟩$444three body interactions
$⟨{\Psi }_{IJK}^{ABC}|$	$\text{contraction}\left[1.0ex\right]{\hat{W}}_N{\hat{T}}_1\text{contraction}\left[1.5ex\right]{\hat{W}}_N{\hat{T}}_1{\hat{T}}_2\text{contraction}\left[2.0ex\right]{\hat{W}}_N{\hat{T}}_1{\hat{T}}_2$ ${\hat{W}}_N{\hat{T}}_1{\hat{T}}_2{\Psi }_0⟩$444three body interactions

Table 4. Table 4: Collection of optimized Goldstone diagrams of the second iteration correlation energy in a fixpoint iterative post-non-canonical HF CCSD.

Index	Contraction	Diagram
A	$\text{contraction}\left[1.0ex\right]{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N\text{contraction}\left[1.5ex\right]{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N\text{contraction}\left[2.0ex\right]{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N\text{contraction}\left[2.5ex\right]{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N{\hat{F}}_N$
B	$\text{contraction}\left[1.0ex\right]{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[1.5ex\right]{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[2.0ex\right]{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[2.5ex\right]{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N$
C	$\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[2.0ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[2.5ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N$
D	$\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[3.0ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N\text{contraction}\left[3.5ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N\text{contraction}\left[2.0ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N\text{contraction}\left[2.5ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N$
E	$\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[2.0ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[2.5ex\right]{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{V}}_N$
F	$\text{contraction}\left[1.0ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[1.5ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N\text{contraction}\left[2.0ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N\text{contraction}\left[2.5ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[3.0ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N\text{contraction}\left[3.5ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[4.0ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N\text{contraction}\left[4.5ex\right]{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N{\hat{V}}_N{\hat{V}}_N{\hat{F}}_N{\hat{F}}_N{\hat{V}}_N$

Equations

$$\hat{T}={\hat{T}}_1+{\hat{T}}_2+\dots ,$$

$${\hat{T}}_{\nu }=\frac{1}{(\nu !{)}^2}\sum _{\begin{array}{c}{i}_1\dots i_{\nu }\\ a_1\dots a_{\nu }\end{array}}{t}_{i_1\dots i_{\nu }}^{a_1\dots a_{\nu }}{\hat{a}}_{a_1}^{†}\dots {\hat{a}}_{a_{\nu }}^{†}{\hat{a}}_{i_{\nu }}\dots {\hat{a}}_{i_1}.$$

$$i,j,k,\dots$$

$$a,b,c,\dots$$

$$p,q,r,\dots$$

$${\hat{H}}_N$$

$$=\sum _{pq}{f}_q^p{\hat{a}}_p^{†}{\hat{a}}_q+\frac{1}{4}\sum _{pqrs}{v}_{rs}^{pq}{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r,$$

$$f_q^p$$

$$v_{rs}^{pq}$$

$$t_{pq\dots }^{rs\dots }$$

$$v_{pq}^{rs}$$

$$⟨{\Psi }_0⟩e^{-\hat{T}}{\hat{H}}_N{e}^{\hat{T}}{\Psi }_0$$

$$⟨{\Psi }_{I\dots }^{A\dots }⟩e^{-\hat{T}}{\hat{H}}_N{e}^{\hat{T}}{\Psi }_0$$

$$\text{contraction}[1ex]E_{\text{corr}}⟵⟨{\Psi }_0|\hat{V}{}_N\hat{T}\text{contraction}[1.5ex]E_{\text{corr}}⟵⟨{\Psi }_0|\hat{V}{}_N\hat{T}\text{contraction}[2ex]E_{\text{corr}}⟵⟨{\Psi }_0|\hat{V}{}_N\hat{T}\text{contraction}[2.5ex]E_{\text{corr}}⟵⟨{\Psi }_0|\hat{V}{}_N\hat{T}{E}_{\text{corr}}⟵⟨{\Psi }_0|{\hat{V}}_N{\hat{T}}_2{\Psi }_0⟩,$$

$$\begin{gathered} \text{contraction}[1ex]⟨{\Psi }_0|\hat{V}{}_N\hat{T}\text{contraction}[1.5ex]⟨{\Psi }_0|\hat{V}{}_N\hat{T}\text{contraction}[2ex]⟨{\Psi }_0|\hat{V}{}_N\hat{T}\text{contraction}[2.5ex]⟨{\Psi }_0|\hat{V}{}_N\hat{T}⟨{\Psi }_0|{\hat{V}}_N{\hat{T}}_2{\Psi }_0⟩=\frac{1}{16}\sum _{pqrs}\sum _{ijab}{v}_{rs}^{pq}{t}_{ij}^{ab}\cdot \\ [\text{contraction}[1.0ex]\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s\hat{a}{}_r\}\{\hat{a}\text{contraction}[1.5ex]\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}\hat{a}{}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}\hat{a}\text{contraction}[2.0ex]\{{\hat{a}}_p^{†}\hat{a}{}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}\hat{a}\text{contraction}[2.5ex]\{\hat{a}{}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}{\hat{a}}_j\hat{a}\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}{\hat{a}}_j{\hat{a}}_i\}+\text{contraction}[1.0ex]\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s\hat{a}{}_r\}\{{\hat{a}}_a^{†}\hat{a}\text{contraction}[1.5ex]\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}\hat{a}{}_s{\hat{a}}_r\}\{\hat{a}\text{contraction}[2.0ex]\{{\hat{a}}_p^{†}\hat{a}{}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}\hat{a}\text{contraction}[2.5ex]\{\hat{a}{}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}{\hat{a}}_j\hat{a}\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}{\hat{a}}_j{\hat{a}}_i\}+ \\ \text{contraction}[1.0ex]\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s\hat{a}{}_r\}\{\hat{a}\text{contraction}[1.5ex]\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}\hat{a}{}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}\hat{a}\text{contraction}[2.0ex]\{{\hat{a}}_p^{†}\hat{a}{}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}{\hat{a}}_j\hat{a}\text{contraction}[2.5ex]\{\hat{a}{}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}\hat{a}\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}{\hat{a}}_j{\hat{a}}_i\}+\text{contraction}[1.0ex]\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s\hat{a}{}_r\}\{{\hat{a}}_a^{†}\hat{a}\text{contraction}[1.5ex]\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}\hat{a}{}_s{\hat{a}}_r\}\{\hat{a}\text{contraction}[2.0ex]\{{\hat{a}}_p^{†}\hat{a}{}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}{\hat{a}}_j\hat{a}\text{contraction}[2.5ex]\{\hat{a}{}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}\hat{a}\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s{\hat{a}}_r\}\{{\hat{a}}_a^{†}{\hat{a}}_b^{†}{\hat{a}}_j{\hat{a}}_i\}] \\ =\frac{1}{16}\sum _{ijab}\left[v_{ab}^{ij}-v_{ba}^{ij}-v_{ab}^{ji}+v_{ba}^{ji}\right]t_{ij}^{ab}=\frac{1}{4}\sum _{ijab}{v}_{ab}^{ij}{t}_{ij}^{ab}. \end{gathered}$$

$${\forall }_{i\in \mathbb{O}}$$

$${\forall }_{a\in \mathbb{V}}$$

$$\text{contraction}\hat{a}{}_i^{†}\hat{a}{\hat{a}}_i^{†}{\hat{a}}_j$$

$$\text{contraction}\hat{a}{}_a\hat{a}{\hat{a}}_a{\hat{a}}_b^{†}$$

$$\underset{\vec{g}(\vec{x})=\vec{0}}{\min }f(\vec{x})\text{ with }f:{\mathbb{G}}^n\to \mathbb{R}$$

$$19^{14}\cdot 6^3\cdot 2^3\approx 1.38\cdot 10^{21}.$$

$$t_i^a=t_{ij}^{ab}=0{\forall }_{ijab}.$$

$$t_I^a$$

$$t_{IJ}^{Aa}$$

$$⟨{\Psi }_I^A⟩\bar{H}{\Psi }_0$$

$$⟨{\Psi }_{IJ}^{AB}⟩\bar{H}{\Psi }_0$$

$$f_I^A=0{\forall }_{I\in \mathbb{O},A\in \mathbb{V}}\text{ with }\mathbb{O}\cap \mathbb{V}=\varnothing$$

Full text

I Introduction

Coupled-Cluster (CC) theory was originally formulated in the field of nuclear physics by Coester and Kümmel Coester_1958 ; Coester_2_1958 . The development of diagrammatic representations of the underlying terms was pioneered by Čížek and Paldus Cizek_1966 ; Paldus_1975 ; Paldus_1977 . These representations were based on the original diagrammatic formulations of Goldstone Goldstone_1957 or Hugenholtz Hugenholtz_1957 , who generalized Feynman’s diagrammatic representation of particle interactions Feynman_1949 , for the many-body perturbation series.

The efficient generation of Hugenholtz- or Goldstone-type diagrams is an ongoing task of computational quantum chemistry and physics, where plenty of algorithms have been proposed in the literature (see e.g. Paldus_1973 ; Csepes_1988 ; Lyons_1994 ; Derevianko_2002 ; Mathar_2007 ). Due to tremendous difficulties in the derivation of higher order ($5,6,\ldots$) CC equations, using algebraic terms and Wick’s theorems Wick_1950 , fast equation derivation is one of the key applications of diagrammatic CC techniques Kallay_2001 , where a lot of progress has been made in the last 20 years (see e.g. Harris_1999 ; Kallay_2001 ; Dzuba_2009 ).

With the advancing algorithm progress deriving CC equations from diagrammatic representations, more and more complex structures in storing and processing the latter representations were created, while the original intuitively understandable diagrammatic structures were lost. These structures however, could be of great importance in the fundamental CC research by supplying insight into more complex contraction patterns or by enabling a sophisticated analyzing tool for new CC approaches.

Considering CC diagrams of higher orders, a large variety of possibilities to correctly illustrate (route) all diagrammatic fragments exists. In particular diagrams, which are not present in the original quantum chemistry CC theory (e.g. including powers of the Hamiltonian ($\hat{H}_{N}^{2},\ldots$) or three particle interrelating operators ($\hat{W}_{N}$)), possess neither unique nor well defined optimal layouts. Finding these can be a very complicated task when done by hand.

In this paper, an automatic Goldstone diagram optimization algorithm („AutoDiag“) is presented. In the current version, it is capable of transforming arbitrary Goldstone diagram topologies to an optimized, readable diagram layout (its representation) meeting a set of (hard) physical and (soft) optical constraints.

As explicitly illustrated in figure 1, diagram topologies are generated from operator expressions using the term generation engine “sqdiag” sqdiag1 ; sqdiag2 while diagrammatic layouts are plotted using the text driven diagram assembling tool “Diag2PS”. Diag2PS

II Theory

In this section, the theoretical basis for CC and its representation using diagrammatic techniques is briefly introduced (for further information see e.g. Shavitt_Book ; Harris_Book ; Crawford_Book ).

II.1 Many body operators in second quantized form

In the CC ansatz, the exact many-particle wavefunction $\ket{\Psi}$ is obtained via the application of the wave operator $e^{\hat{T}}$ to the Hartree-Fock (HF) reference determinant $\ket{\Psi_{0}}$. Here, $\hat{T}$ denotes the total cluster operator with

[TABLE]

where the individual Tν operators denote $\nu$-particle substitutions via

[TABLE]

Using the common index notations for occupied ($\mathbb{O}$) and virtual ($\mathbb{V}$) orbital spaces

[TABLE]

the normal ordered Hamiltonian $\hat{H}_{N}$ takes the form

[TABLE]

where $f_{q}^{p}$ and $v^{pq}_{rs}$ denote the one- and two-particle integrals

[TABLE]

respectively. In general, the amplitudes $t$ and the integrals $v$ are antisymmetric with respect to index permutations such that

[TABLE]

II.2 Algebraic evaluation

The projected and similarity transformed CC equations are given by

[TABLE]

where capital indices $A\ldots$ and $I\ldots$ denote external projections.

The similarity transformed Hamiltonian $\bar{H}=e^{-\hat{T}}\hat{H}_{N}e^{\hat{T}}$ can be expanded in the Baker-Campbell-Hausdorff (BCH) formula Baker_1902 ; Campbell_1897 ; Hausdorff_1906 to yield nested commutators of $\hat{H}_{N}$ and $\hat{T}$, which naturally vanish for quintuply (and higher) nested commutators due to the two-particle interrelating nature of $\hat{H}_{N}$. Inserting the BCH series of $\bar{H}$ into the CC equations (5) and (6), only fully contracted Wick-terms survive. Furthermore, all participating operators from $\hat{H}_{N}$ must be at least singly contracted to all $\hat{T}$ operators on its right.

As an example for algebraic evaluation consider the term

[TABLE]

which contributes to the correlation energy (5). To algebraically evaluate (7), all possible full contraction patterns of the underlying second quantized operator strings must be collected. These result in products of Kronecker-deltas, which bring down the general index summation ($p,q,r,s$) by binding them to particle/hole indices ($i,j,a,b$) of the amplitudes $t$. Using antisymmetry, all redundant two-particle integrals $v$ can then be merged to end up with

[TABLE]

As is clearly recognized, even this fairly simple algebraic derivation example turns out to be annoyingly complex. In part, this complexity arises from the fact that the Wick theorems do not incorporate the antisymmetry of the tensors $v$ and thereby produce a lot of redundant terms.

II.3 Diagrammatic evaluation

A different approach, pioneered by Čížek and Paldus Cizek_1966 ; Paldus_1975 ; Paldus_1977 , was found when using diagrammatic techniques to derive the CC equations. In the particle/hole formalism (derived and applied also in the algebraic framework), the definition of all second quantized operators acting on $\mathbb{O}$ is inverted to yield the „quasi-particle operators“ $b$:

[TABLE]

These may now be illustrated by directed lines in a Feynman-type space-time diagram (c.f. fig. 2).

Here, hole-operators (acting on $\mathbb{O}$) are represented as directed lines facing down (moving backwards in time) and particle-operators (acting on $\mathbb{V}$) are represented as directed lines facing up (moving with time). A further distinction is made for creators and annihilators with respect to a certain event in time (e.g. a particle interrelating operator contracted to the individual particle/hole lines). All outgoing particle/hole lines represent creators ($\hat{a}_{i}^{\dagger}$ and $\hat{a}_{a}^{\dagger}$) and all incoming particle/hole lines represent annihilators ($\hat{a}_{i}$ and $\hat{a}_{a}$).

The most important feature of the particle and hole line representation is found when considering Wick-contractions. The only surviving contraction types are

[TABLE]

As it turns out, these contraction rules are perfectly satisfied when connecting particle/hole lines heads to tails (c.f. fig. 3). First of all, it is only possible to connect two arrows if both of them are facing up or facing down allowing for contractions among two holes or two particles only. Secondly, an outgoing line must always be connected with an incoming line thus only ensuring contractions among creators and annihilators.

To fully represent particular many-body integrals of the CC equations (5) and (6) – as e.g. illustrated in figure 4 – the occuring particle interrelating operators of the Hamiltonian $\hat{H}_{N}$ and the cluster operator $\hat{T}$ need to be represented. In general, particle interrelating operators are illustrated by horizontally arranged vertices connected by a dashed interaction line, which are connecting several particle/hole lines.

The individual cluster operators $\hat{T}_{\nu}$ have fixed particle/hole indices and always promote particles from $\mathbb{O}$ to $\mathbb{V}$. Therefore, they always possess the same „V“-shape of connected particle/hole lines. Cluster operators $\hat{T}_{\nu}$ with $\nu\geq 2$ are illustrated by $\nu$ horizontally arranged vertices connected by a bold cluster line.

In figure 4, the same example (eq. 7) as evaluated in section II.2 is illustrated using the corresponding Goldstone representations. From the definition of the integral, the participating operator representations (with their correct ranks) may be extracted and illustrated according to their time of operator application (with respect to the ket Fermi-vaccuum $\ket{\Psi_{0}}$) in a vertically advancing time scale. Connecting all particle/hole lines as mentioned above (heads to tails) leads to the complete Goldstone diagram, which consists of two particle/hole line loops in this example. This diagram is completely translatable (including sign and prefactor) to the algebraic evaluation of the integral.

III Discrete optimization problem: Diagram routing

Since there is no hard criterion for an optimal Goldstone diagram layout, a manageable global optimization problem of type

[TABLE]

is defined. Here, $\vec{x}$ represents the target quantity, wich completely defines one specific layout (e.g. all information of solution space $\mathbb{G}$ to fully characterize all $n$ vertices).

In order to guarantee physically correct Goldstone diagrams, i.e. a conserved operator time ordering, specific constraints $\vec{g}(\vec{x})$ are placed upon the minimalization problem (c.f. subsection III.3). This leads to a significantly reduced area of $\bar{\mathbb{G}}^{n}\subseteq\mathbb{G}^{n}$ where all constraints are fulfilled. This remaining area represents the solution space the algorithm should traverse.

Finally, for each layout $\vec{x}$ a certain cost value $f(\vec{x})$ (c.f. subsection III.5) is defined, such that the global minimum of $f$ represents an optimally routed diagram.

III.1 Explicit spatial discretization

The automatic routing of Goldstone diagrams clearly resembles a high-dimensional problem, since all vertices of the latter must be placed in their optimal position. Therefore, two dimensional diagrams without any constraints possess a dimensionality of $2n$ (for $n$ vertices with $x$ and $y$ coordinates). Per definition however, the optimization problem is not discrete. Only when introducing a discrete solution space for each diagram, an optimization algorithm can perform well. This work employs the text driven diagram assembling tool „Diag2PS“ Diag2PS , which is capable of illustrating Goldstone, Hugenholtz or classical Feynman diagrams. In this tool, a discrete grid of thirds is used to place the vertices of the assembled diagrams. This grid was used to discretize the vertex positions in this work.

An example illustration of the grid-based Goldstone diagram optimization problem is shown in figure 5. The diagram on the left hand side was generated randomly following its time-ordered topology from the CCSD expansion, to illustrate a random starting point of the optimization. On the right hand side, the desired minimum is shown. Each vertex position is exactly characterized by its integer $x$ and $y$ coordinates on the grid.

III.2 PHL structure

Knowledge about the coordinates $x$ and $y$ alone is not sufficient to fully characterize Goldstone diagrams. It is also mandatory to optimize the topology, i.e. the particle/hole line (PHL) structure, of the diagram. This structure is composed in the following fashion:

• (i)

A single PHL is defined as a directed line connecting two distinct vertices. The spatial extent (curvature, straightness) of the latter is arbitrary as long as the direction (facing up or facing down) as well as the starting and ending vertices of the PHL are conserved.

• (ii)

Several PHLs may form a connected segment (a path) of alternating incoming and outgoing PHLs advancing over several vertices (c.f. fig 6). Therefore, every PHL segment can either start and end at exactly the same vertex (representing a loop, c.f. segments (a) and (b)) or start and end at different external vertices (c.f. segment (c)), which stop the connected path.

• (iii)

The collection of all PHL segments of a Goldstone diagram is denoted by its PHL structure, which fully characterizes its topology.

The PHL structure of any Goldstone diagram is naturally discrete since (i) all PHLs must start and end at specific single vertices and (ii) all incoming or outgoing lines are uniquely mapped to exactly one vertex (i.e. no vertex can be connected by more than one incoming or outgoing line).

III.3 Constraints

In case of Goldstone diagrams there are several constraints for a physically correct representation of the latter. In summary, the following constraints are placed on the optimization problem:

• (i)

Every interaction or cluster vertex/line (representing one specific particle interrelating operator from the Hamiltonian or one specific cluster operator) is identified with one distinct relative position in time. This time scale is based on the order of operator applications to the ket Fermi vaccuum $\ket{\Psi_{0}}$, where first acting operators are illustrated earlier in time then last acting operators. In this work, the time axis in Goldstone diagrams is illustrated as advancing from the bottom to the top.

• (ii)

Since interaction or cluster lines can consist of several vertices (and each interaction or cluster line is only identified with one time position), all those vertices must possess the same time position and must therefore be illustrated horizontally arranged.

• (iii)

Whenever two operators commute (e.g. several $\hat{T}$ operators for disjoint $\mathbb{O}$ and $\mathbb{V}$) they are capable to interchange their time position arbitrarily. Conventionally, such operators are illustrated at the same time position to illustrate their interchangeability.

• (iv)

An exception is made for all external lines originating from projections. In principle, external projections always possess the highest position in virtual time since they are always acting from the left ($\bra{\Psi_{I\ldots}^{A\ldots}}\hat{X}=\bra{\Psi_{0}}\hat{a}_{I}^{\dagger}\ldots\hat{a}_{A}\hat{X}$). However, since they are represented by isolated particle/hole lines, they are allowed to fall down to the individual cluster or interaction lines they are connected to. As a hard constraint they must never fall below their connected interaction or cluster lines.

Taking all of these restrictions for granted, the resulting Goldstone diagrams can always be correctly translated to their algebraic form. With increasing diagram size however, the number of physically correct representations increases enormously while only few are still intuitively readable (analyzable). This does not only include the spatial positioning of all vertices but also the PHL structure.

III.4 Degrees of freedom

If all constraints are fulfilled, the following degrees of freedom allow for an optimization of the cost:

• (i)

The horizontal axis of a Goldstone diagram has no algebraic implications. If a vertex is illustrated more to the left or to the right for instance, has no effect on the algebraic correctness of the diagram.

• (ii)

The actual position in time of any Goldstone interaction or cluster line is arbitrary as long as it does not violate the time ordering itself. This means that vertices may be placed closer together or further apart in the vertical time dimension without interchanging their position.

• (iii)

Due to the antisymmetry of the included tensors (c.f. equations 3 and 4), pairs of two outgoing (creators) or two incoming (annihilators) PHLs connected to the same interaction or cluster line may be interchanged. This changes the topology of the Goldstone diagram and may therefore lead to a significant change in the PHL structure, which e.g. includes the number of loops of a diagram (c.f. fig. 7).

III.5 Cost evaluation

The key concept for the global optimization problem is the definition of a cost function $f$ evaluating the quality of individual diagrams. In general, low costs should correspond to a good diagrammatic quality and high costs to a diagram of poor quality. For Goldstone diagrams, several drawing conventions were defined, which are fulfilled by diagrams that are considered more elegant than others. These conventions should, in the first place, increase the readability of the individual diagrams and should somehow imitate human-drawn diagrams (i.e. should prefer symmetric diagrams or prefer diagrams without line crossings, etc.).

III.5.1 Static costs

Whenever evaluating Goldstone diagrams, there are static cost factors appearing solely through the individual fragments of the diagram. Those static costs are summarized in table 1 for different example diagrams. The individual cost values for each of the illustrated factors were chosen empirically.

1. (A)

Since any form of loop is easily recognizable and is beneficial for a scaling analysis, the algorithm should always try to find as many loops as possible for the given diagram. Every loop therefore lowers the total costs by $10$ units. 2. (B)

Crossings of any kind on the other hand can be misleading and should be avoided. This is ensured by adding $5$ units per PHL crossing to the total costs. 3. (C)

Since interaction or cluster line crossings can only occur horizontally, a crossing of e.g. two $\hat{T}_{2}$ cluster lines may easily be misinterpreted by a $\hat{T}_{4}$ cluster line. Therefore, any interaction or cluster line crossing is rated worse than any PHL crossing by applying twice the cost of a PHL crossing. 4. (D)

If two lines are close to each other, it is harder to distinguish between them. In such a case, $1.5$ units are added to the total costs forcing all lines to lie within a certain minimal distance. 5. (E)

Furthermore, $100$ units are added whenever two vertices are arranged at the same position, which can also lead to irritating diagrams. In the illustrated case, two external vertices fall together forming a joined vertex, which could be wrongly recognized as a $\hat{T}_{1}^{\dagger}$ vertex. 6. (F)

Finally, a symmetry check is executed evaluating the existence of a plane of reflection in the center of the diagram. If such a plane was found, one unit is subtracted from the total costs.

III.5.2 Dynamic costs

In addition to fixed costs, several dynamic (variable) costs are included in the implemented cost function.

• (I)

To ensure a reasonable line length, all particle/hole and interaction line lengths are calculated during the cost evaluation. For straight PHLs, a preferred line length of $\sqrt{\frac{10}{9}}$ was defined, which corresponds to the length of a line advancing one full unit (three thirds in the solution grid) upwards and one third of a unit sideways. For PHLs, which are part of a loop, the preferred line length corresponds to the difference in time-orders of the two interaction lines connected by the loop. Any absolute deviation of the calculated line lengths to the preferred ones is added to the total costs.

• (II)

Whenever two particle or two hole lines are connected to the same vertex, a straightness criterion is checked. If the found particle or hole lines are not collinear, the deviation in $x$-direction to collinearity is added to the total costs.

IV Approach of solution: Genetic Algorithm

The Goldstone diagram routing problem as defined in section III as a global minimization of the cost function $f$ represents a discrete high-dimensional optimization problem with exponentially increasing complexity with the number of vertices. This suggests the application of a genetic algorithm as outlined within the next subsections.

IV.1 General structure

First introduced by Holland Holland_1975 , genetic algorithms resemble a particular field of algorithms that use techniques inspired by evolutionary biology. These techniques include selection, crossover (and therefore inheritance), mutation and death. Genetic algorithms are considered as efficient solution algorithms for high-dimensional, global optimization problems in a discrete solution space, Kumar_2010 which have been used in CC related research fields (see e.g. Engels_2011 used by Hanrath_2010 ).

Figure 8 shows a sketch of the steps involved a genetic algorithm as pseudo code.

Starting with the creation of random individuals, an initial guess of $n$ individuals is instantiated. In terms of genetics, these individuals are recognized as genes, that make up the first population. Each individual gene is rated according to its quality with respect to the problem in question. Genes of higher quality may be rated with a lower cost (or a higher fitness) than genes of poorer quality. Certain genes may then be selected as suitable mates, creating a new gene, which can be generated as a crossover of the mating genes. This step of creating genetic children is called crossover. To ensure genetic diversity, every genetic child has a specific chance (or rate) to undergo a mutation process, resembling a random change of its contents (its genetic material). Via these crossover and mutation processes, new genes are derived from already existing ones while keeping a controlled degree of random variation through the employed mutation rate. Especially for global optimization problems, a key concept of solution algorithms is the ability to converge across local minima into the desired global minimum. Even if the local minima are well spread within the whole solution space, the employed mutation process may grant this ability depending on the number of generations, the population size and the mutation rate.

After introducing all generated children into the population, the latter undergoes a selection process, where the weakest genes are erased from it. This keeps the population size fixed to $n$ genes and resembles an important optimization step, since all inaccurate genes are removed and are therefore not allowed to pass their genetic material to possible children. Finally, if no further evolutionary progress was made (e.g. no improvement was found) the algorithm terminates with the final result being found in the best gene of the final population. In every other case, the population will undergo all steps again, starting with the crossover of genes.

Although the algorithm itself has a rather fixed structure, the exact formulation of crossover and mutation are strongly dependent on the underlying optimization problem and are crucial for the performance of the algorithm.

IV.2 Reproduction

The reproductive process (crossover and mutation) resembles one of the most crucial steps of any genetic algorithm towards convergence. In particular for Goldstone diagrams, reproduction is neither trivial nor in some fashion canonically defined. As outlined in the next subsections, a genetic crossover (c.f. subsection IV.2.1) based on diagram fragmentation was found to perform well while different mutational processes sorted into two mutational phases are employed (c.f. subsection IV.2.2).

IV.2.1 Crossover

A definition of a reasonable genetic crossover of two Goldstone diagrams was found employing a fragment-based ansatz. The principle concept of this ansatz is illustrated in figure 9.

In this ansatz, two parent Goldstone diagrams produce a child diagram via fragmentation, competition and reassignment. The individual parent diagrams are illustrated at the top of figure 9. Both parents are fragmented into PHL segments

(c.f. subsection III.2) by breaking the participating interaction or cluster lines. Those fragments are represented in solid and dashed frames, illustrating their parent affiliation. All corresponding fragments are then competing against each other, where the costs of the individual fragments are evaluated and compared.

For this fragmentation and cost evaluation, the parametrization of PHLs into individual segments is beneficial. Since all PHL segments are fragments per definition, no explicit fragmentation process is needed. Iterating over all vertices contained in one PHL segment allows for an efficient fragment cost calculation, which is needed to compare the individual fragments later on. The fragment cost calculation itself is not different from a usual cost calculation as described in subsection III.5 checking all the proposed criteria only for the selected subset of all vertices (with a few obvious exceptions e.g. symmetry or length costs for interaction lines).

In the final reproductive step, the best fragments are reassembled to result in the final child diagram. This is done by directly copying the $x$ and $y$ coordinates of the winning fragments to the child.

There are two major issues concerning the fragment based crossover ansatz:

1. (i)

The fragment based crossover is not always possible. For any successful reassignment, the fragments of both parents must be exactly equal to each other. Using the antisymmetry principle, the outgoing or incoming lines of any interaction or cluster line may be interchanged, thus creating different fragments. In any such cases, the two parent costs are directly compared to each other and the better parent diagram will be cloned. 2. (ii)

The fragment’s relative coordinates towards the other fragments are not taken into account. Therefore, many crossovers may lead to poor diagrams with equal vertices or many line crossings. With increasing generation however, only diagrams that are beneficially reproducible survive the algorithm. In fact, a direct comparison of the simple reassignment to a more sophisticated reassignment procedure revealed that much better results and a better efficiency were obtained when using the proposed simple approach. The collective intelligence of the genetic algorithm seems to reside in the selection process and not within a sophisticated crossover for this particular problem.

IV.2.2 Mutation

The mutation process assures the needed genetic diversity to inspect all possible areas of the solution space. It is responsible for a successful convergence of the algorithm to the global minimum. To ensure this diversity, the mutation process must be capable of reaching the global minimum by accident within a single mutation of any diagram. In particular, this includes the topology altering usage of antisymmetry (c.f. subsection III.4) throughout the whole algorithm.

Additionally to an exponentially decreasing mutation rate, it was found benefical (faster converging) to define different mutational phases I and II, where different operations are carried out:

Phase I: Large diversity via random absolute positioning

Phase I features mutations of large diversity, where $x$- and $y$ coordinates of all vertices are randomly assigned to new values with the chance of the exponentially decreasing mutation rate.

Phase II: Small diversity via differential updates

Phase II features minor mutations, which employ differential updates to the vertex positioning. All possible mutations during phase II are illustrated in table 2.

1. (A)

With the given mutation rate, certain PHL fragments are selected for a random movement in $x$ direction between $-1$ and $+1$ units. All other vertices of all other fragments are kept constant. This may cure too long or too short interaction or cluster lines. 2. (B)

A scissoring movement was employed to account for large interaction or cluster lines, where all vertices are not equally distributed. Unique pairs of PHL segments are chosen with the chance of the given mutation rate. These pairs are then moved towards or away from each other with $\{-1,\ldots,+1\}$  (randomly chosen) units in $x$-direction. 3. (C)

An interaction line shift was employed, moving all vertices associated with randomly chosen interaction or cluster lines by $\{-2,\ldots,+2\}$  (randomly chosen) units in $x$-direction and by $\{-1,\dots,+1\}$ units (if time order conserving) in $y$-direction. 4. (D)

Finally, a flip of individual PHL segments was employed, were all vertices of a certain segment are mirrored with respect to a central plane of reflection. This mutation is particularly useful for diagrams containing $\hat{F}_{N}$ vertices, where the cross is oriented in an unfavorable direction.

V Convergence analysis

In the following subsections, several analyses of the proposed algorithm are demonstrated.

Subsection V.1 presents an an in-depth view on convergence and complexity of the genetic algorithm. An arbitrary example of large complexity from the fixpoint iterative CCSD expansion is used.

The reproductive process is analyzed in subsection V.2 while timings of the algorithm are given in subsection V.3.

V.1 Complexity and convergence behavior

In figure 10, the optimization route explicitly showing the leading diagrams throughout the optimization, is illustrated for the optimization of the diagram constructed in figure 16.

This particular optimization was conducted using the default settings of a population size of 40 genes times the number of vertices, which corresponds to $p=560$ in the illustrated case. The convergence threshold of the algorithm was set to a cost difference from worst to best diagram of $10^{-6}$ units. The mutation rate was chosen to exponentially decay from $1.0$ to $0.1$ within 5600 generations ($g=10p$) (if not converged earlier). In case of no convergence after the 5600 generations the algorithm propagates with a constant mutation rate of $0.1$. As described in subsection IV.2.2, two mutational phases were employed to reduce the number of generations needed to converge. Phase I (larger mutations) was applied from generation 1 to 2799 and phase II (weaker mutations) was applied from generation 2800 ($\frac{g}{2}$) to convergence.

Figure 10 (a) shows the best diagram of the initial population (created randomly) and (f) the converged diagram. Clearly, the initial diagram at generation 1 represents a poor solution guess with many PHL crossings. However, due to the used parametrization (c.f. section III), all diagrams automatically follow the correct virtual time ordering. This is even true for the first generation only consisting of randomly created diagrams. Within the large mutation range from generation 1 to 2799 to diagram (b), the algorithm proceeds in avoiding the fundamental costs such as e.g. crossings or line vicinities. In the small mutation range afterwards through diagrams (c) to (f), the optimal routing is constructed, applying only small shifts of single vertices or interaction lines.

The optimized diagram contains 14 vertices with degrees of freedom in their $x$-coordinates and three interaction line collections with unique time positions, which possess degrees of freedom in their $y$-coordinates. Therefore, the optimization problem contains $17$ positional degrees of freedom. The PHL structure for each $\hat{V}_{N}$ contributes two ($(2!)^{2}-2!=2$) additional degrees of freedom.

The minimal solution grid of thirds contains $19\times 6$ unique vertex positions. Since there is no restriction that different vertices can not occupy the exact same position, the size of the possible solution space (neglecting time ordering) is exceedingly large with roughly

[TABLE]

Due to this large solution space, any brute force algorithm is impossible and the solution via a genetic algorithm in only 4558 iterations (generations) presents an efficient solution.

To further analyze the convergence behavior, the optimization course was illustrated as a function of the individual diagram costs in figure 11.

In this illustration, every gene of every population (in every generation) is described by a single point corresponding to its costs. Therefore, it is in principle possible to track how individual genes propagate through the genetic algorithm from their first appearance to their potential death.

During phase I, a fast decrease in costs within the first generations is recognizable. Subsequently, the costs decrease linearly up to the end of phase (I). The cost range from best to worst diagram is kept nearly constant. Therefore, the genetic diversity of all populations is large enough to ensure that the algorithm does not converge too early in any local minima.

In phase II however, the cost range is drastically reduced. The removal of phase II lead to exceedingly large generation numbers needed to converge the algorithm. As a compromise, the phase transition was implemented at $g=10n/2$, which produced optimal results for every tested diagram.

V.2 Crossover and mutation ratio

Children are created by directly combining the coordinates of the winning fragments from their parents. This may lead to diagrams which are worse compared to their parents. As stated in subsection IV.2.1, the fragment-based crossover for Goldstone diagrams is not always capable of producing a genetically altered child. In any such cases, the algorithm will proceed by creating a clone of one of its parent diagrams.

Using the same example of the diagram constructed in figure 10, the crossover and mutation processes were analyzed. This was done by calculating the number of diagrams improved by crossover (with respect to both their parents) and by mutation, purely by crossover, purely by mutation and by none of both. The percental ratios from all created children during each generation are illustrated in figure 12.

Clearly, pure crossover ($20\%$) dominates over combined crossover and mutation and pure mutation during the beginning of the algorithm. This effect results from the large mutation rate and range at the beginning of the algorithm. While mutation keeps a large genetic diversity in the algorithm (by producing diagrams in every area of the solution space), pure crossover can produce an improved child in nearly every fifth case. Towards the end of phase I, pure crossover slowly decreases below 10% until phase II is entered.

Due to the smaller mutation range in phase II, improvement by pure crossover starts increasing instantly while improvement by pure mutation enters with roughly 5% to 10% of improved children. While the improvement by pure mutation stays low around the 5% mark, the improvement by pure crossover decreases again until the final converged diagram is constructed at generation 4558. After this point pure crossover and pure mutation tend towards 0% since there is no additional cost gain from altering the latter diagram.

A similar picture may be obtained for different diagrams of comparable complexity and it is therefore possible to conclude that the fragment-based crossover works.

V.3 Timings

The computing time of the genetic algorithm was statistically analyzed for a large collection of diagram optimizations. These include different CC approaches including different Hamiltonians (normal, three-particle interrelating and $\hat{H}_{N}^{2}$) as well as variational and similarity transformed CC approaches. The time ordered topologies were obtained by the term generation engine “sqdiag” sqdiag1 ; sqdiag2 and the optimizations were conducted on a typical office machine using a single core of the Intel Core i3-8100 CPU with 3.60 GHz.

In figure 13, the individual diagrammatic complexities are plotted against the wall clock computing time (on a single computing core) needed to achieve convergence. The complexities were calculated accordingly to subsection V.1.

Clearly, an increased complexity of the solution space in general leads to higher computing times. A large broadening effect is recognizable, where different optimizations of comparable complexity may take different computing times (up to two orders of magnitude). This may result from various starting configurations of the initial population, which may take different numbers of generations to converge.

Due to these different generation numbers, any direct relation of diagrammatic complexities to real computation times is biased. Therefore, the average wall clock time for one generation was calculated and illustrated with respect to the diagrammatic complexity (c.f. figure 14).

In general, wall clock times of under 1 ms to 10 ms per generation were found. Again, a slight broadening effect is visible, which appears in a small range of under 10 ms.

Overall, sufficiently small wall clock times were found such that no conducted optimization took longer than 100 seconds. It is important to state that up to this point no thorough code optimization was conducted. So far the emphasis was on functionality of the algorithm rather than efficiency. Additional code optimizations would presumably lower the total computation time significantly.

VI Application

Several Goldstone diagram optimizations were conducted from different CC approaches to check the functionality of the developed algorithm. This includes example diagrams from

• •

the standard similarity transformed CCSDTQ ansatz (up to quadruple projections),

• •

the variational CC ansatz for CCSDTQ (only linear terms $e^{\hat{T}}\approx 1+\hat{T}$),

• •

the variational CC ansatz for CCSD (linear and quadratic terms $e^{\hat{T}}\approx 1+\hat{T}+\hat{T}^{2}$),

• •

the similarity transformed CCSD ansatz employing $\hat{H}_{N}^{2}$ (up to single projections) and

• •

the similarity transformed CCSDT ansatz for $\hat{H}_{N3}$ (up to triple projections).

The time-ordered diagram topologies were obtained from the term generation engine “sqdiag” sqdiag1 ; sqdiag2 . To illustrate the quality of the obtained diagram layouts, a selected collection is given in the supporting information. All diagram optimizations were conducted multiple times for different initial populations. All of those optimizations, found the same or qualitatively equal solutions even for initial populations of poor quality. To discuss benefits and potential flaws of the developed algorithm, a small collection of diagrams with mixed complexity from all conducted optimizations is illustrated in table 3.

All diagrams were successfully optimized to possess the maximum amount of loops, as few line crossings as possible (particle/hole, interaction and cluster) and no too close PHLs increasing the readability of the latter. Furthermore, more complex loops containing three or more vertices are correctly optimized and illustrated. Even the 4-loop case always containing a PHL crossing (e.g. in the variational CC ansatz) is illustrated with circular PHLs at the outer sides signalizing its loop status. Additionally, symmetric diagrams are preferred (if there is no additional loop by breaking symmetry). We conclude that the developed algorithm is capable of producing essentially optimal Goldstone diagrams.

VI.1 CC fixpoint iterations

To further investigate the functionality of the developed algorithm, the CC fixpoint iterations were recreated in a purely diagrammatic representation. These present a good testing scenario because new diagrams with rapidly increasing complexity are generated only from the standard CC diagrams by diagrammatic substitutions (see e.g. Crawford_Book ). Therefore, no further input from the user is required once all diagram topologies are read in and the program is started. The proposed algorithm is capable of arbitrary Goldstone diagram substitutions (arbitrary CC truncation). Here, the results for CCSD are discussed in detail.

Considering the CCSD equations of the zeroth iteration, all amplitudes are set to zero with

[TABLE]

This of course directly results in a zeroth iteration correlation energy $E_{\text{corr}}^{(0)}$ of zero, since in every CCSD energy diagram either T1 or T2 amplitudes are present:

[TABLE]

The new first order T1 and T2 amplitudes are obtained from the singles and doubles projections (c.f. figure 15). Since all amplitudes are zero, it is

[TABLE]

and therefore the first three terms of the singles and doubles equations as illustrated in figure 15 simplify to

[TABLE]

Algebraic manipulation then leads to the first iteration T1 and T2 amplitudes consisting of the only diagrams not containing any amplitudes (purely external diagrams):

[TABLE]

For post-canonical HF CCSD the first iteration T1 amplitudes are zero due to the diagonalized Fock matrix, for which

[TABLE]

holds.

For the sake of the algorithm testing however, post-non-canonical HF CCSD diagrams were considered to gain an even more rapidly increasing diagram complexity. In figure 16, one particular example of the second iteration CCSD energy expression is constructed.

Table 4 shows a small collection of optimized Goldstone diagrams of different complexity from the second iteration correlation energy. All diagrams occuring in zeroth, first and second iteration (energy, singles and doubles) are illustrated in the supporting information. In perturbation theory, it is customary to keep track of the energy denominators (e.g. $D_{1}$ and $D_{2}$) by the inclusion of horizontal lines in the corresponding Goldstone diagrams. Every PHL crossed by these horizontal lines contributes to the energy denominator. Harris_Book For the particular fixpoint iterative procedure presented here, these lines are not explicitly shown. It is also important to note that the usual rules to translate perturbation theory diagrams to their algebraic expressions (including their energy denominator) do not apply to the diagrams shown here since only those PHLs explicitly part of the substitution contribute to either one $D_{1}$ or one $D_{2}$ denominator. To illustrate the fixpoint substitution, one interaction line (from the original leftmost $\hat{F}_{N}$ or $\hat{V}_{N}$) is drawn at the highest possible relative position in time.

Clearly, all illustrated diagrams (A to F) as well as all fixpoint iterative diagrams in the supportig information are reasonably optimized. Every single one of them is easily readable, analyzable and possesses the physically correct time ordering. The algorithm is capable of the optimization of several loop constructs. These include the „ladder“ 3-loop (c.f. diagram E), the crossed 4-loop (c.f. diagrams A, D, F) as well as the direct 2-loop (c.f. diagrams B, C, D, E, F) or combinations of all. Furthermore, antisymmetry is benefically used to increase the number of loops of all diagrams as e.g. clearly visible for diagram C.

In the inevitable case of a PHL crossing (c.f. diagrams A, D, F) the latter is routed in a still easily recognizable way such that the overall readability is not decreased. Furthermore, symmetric diagrams (if possible) are preferred over non-symmetric ones (c.f. diagram C).

VII Conclusion

Throughout this paper, the concept, analysis and application of the proposed automatic routing algorithm for Goldstone diagrams were illustrated.

Section III focused on the discrete Goldstone diagram optimization problem. A discretization of Goldstone diagrams was achieved in a grid of thirds while PHL structures were defined as a collection of PHL segments (c.f. subsections III.1 and III.2). Suitable constraints for physically correct Goldstone diagrams were established (c.f. subsection III.3) to reduce the possible solution space. The remaining diagrammatic degrees of freedom (c.f. III.4) represent the parameters, that need to be optimized in order to find a well routed global minimum. A cost function (c.f. III.5) evaluating the quality of individual Goldstone diagrams was defined.

Due to the exponentially increasing complexity of the proposed optimization problem, a genetic algorithm was developed. A successful crossover approach was implemented in form of a fragment-based crossover (c.f. IV.2.1) transferring genetic information from two parent diagrams to a newly created child diagram. Different mutation processes (c.f. IV.2.2) were implemented to keep genetic diversity throughout the algorithm and prevent premature convergence effects.

In section V, the convergence behavior, including a comparison of diagram complexities and the implemented crossover and mutation processes were analyzed. It was shown, that Goldstone diagram solution spaces can possess enormously large complexities (c.f. V.1) paying justice to the developed genetic algorithm to reach convergence within a reasonable amount of iterations. Via the investigation of crossover and mutation processes (c.f. V.2), it was possible to show that the fragment-based crossover is beneficial by creating improved child diagrams compared to their individual parent diagrams. Improvements by mutation processes alone primarily occur during the final generations of the Goldstone routing algorithm.

Section VI dealt with the applicability of the developed algorithm to different CC approaches. It was possible to reach diagrams, which fulfill all criteria described in subsection III.5, in every conducted optimization. As a more complex testing ground, the fixpoint iterative substitution of Goldstone diagrams for the standard similarity transformed CCSD in the non-canonical post-HF procedure was evaluated in subsection VI.1. Additional applications of the proposed algorithm may be found in the supporting information.