Parallelized traveling cluster approximation to study numerically spin-fermion models on large lattices

Lattice spin-fermion models are important to study correlated systems where quantum dynamics allows for a separation between slow and fast degrees of freedom. The fast degrees of freedom are treated quantum mechanically while the slow variables, generically referred to as the “spins,” are treated classically. At present, exact diagonalization coupled with classical Monte Carlo (ED + MC) is extensively used to solve numerically a general class of lattice spin-fermion problems. In this common setup, the classical variables (spins) are treated via the standard MC method while the fermion problem is solved by exact diagonalization. The “traveling cluster approximation” (TCA) is a real space variant of the ED + MC method that allows to solve spin-fermion problems on lattice sizes with up to ${10}^3$ sites. In this publication, we present a novel reorganization of the TCA algorithm in a manner that can be efficiently parallelized. This allows us to solve generic spin-fermion models easily on ${10}^4$ lattice sites and with some effort on ${10}^5$ lattice sites, representing the record lattice sizes studied for this family of models.

Figure 1

Two-dimensional schematic of the TCA approach. Here a finite lattice is displayed with classical DOF at each site, represented by the red arrows. The Hamiltonian defines the coupling of the classical spins with the itinerant electrons on a lattice. The TCA algorithm consists of proposing an update at a particular site (encircled in blue) for the classical spins. The update is accepted or rejected based on the energy of a cluster built around that site, indicated by the shaded (green) rectangle. Here the cluster size is $N_c=9$. In a system sweep, the above procedure is carried out by visiting each site of the system sequentially.

Figure 2

The flowchart for a single sequential Monte Carlo system sweep in the TCA approach. Here, we consider a lattice of $N$ sites. $S\left[\right]$ is an array containing the classical DOF at all the $N$ sites. $C\left[i_s\right]$ is an array of size $N_c$ that reads in the classical DOF from $S\left[\right]$ around the “update site” labeled by $i_s$. The loop over $i_s$ runs over all the $N$ sites of the system, ensuring that the cluster is built around all the sites sequentially and that an update is proposed at each site. The flowchart is discussed in the text.

Figure 3

One-dimensional example using $N=8$ and $N_c=4$. The light-blue sites are the unupdated sites. The red (hatched) site is chosen to be the site where the update is attempted (update site). The clusters are indicated by the boxes. The different rows present the cases of the cluster traveling sequentially from left to right during a single-system sweep. The filled (blue) sites indicate the sites where an update has been attempted. The system and cluster both have periodic boundary conditions. For rows 1 through 5, the clusters are built based on the original (unupdated) classical variable configurations and can be constructed simultaneously instead of serially as in TCA. For rows 6 through 8, the clusters require the results of the update attempts on sites 1 through 3, respectively. So these clusters have to wait till row 4 and then can be built simultaneously.

Figure 4

Parallelized TCA flowchart for a single-system sweep. Unlike Fig. 2, the loop over the total number of sites of the system is split into two loops: the outer one runs over the number of blocks, while the inner one runs from 1 to the ratio of the number of sites in a block ($N_S$) to the number of processors ($N_P$). The gray rectangles, labeled $P{C}_1$, $P{C}_2,...,$ $P{C}_{NP}$, are computed on processors with rank 1 through $N_P$, respectively. The steps in a typical gray block are displayed on the left. See text for discussion.

Figure 5

(a) Time required for asynchronous ED and full ED with message passing against the number of processors $N_P$ at a fixed $U=8.0t$. See text for definitions. The data is presented for a $N=8^3$ lattice using 2000 MC system sweeps at a fixed temperature. The dashed line is the plot of $1/N_P$; (b) Time needed for 200 system sweeps for the full calculation with message passing plotted against $L$, for a system size $N=L^3$. The data for different number of processors, $N_P$, are shown. $N_P$ values are indicated on the right. The dashed line indicates that the computational cost for solving a $N={32}^3$ system using PTCA with 64 processors is almost the same as the time needed using TCA on an $8^3$ system with a single processor. Calculations were done using multiple Intel Xeon E5-2670 processors, which have eight cores with base frequency of 2.6 Ghz and 8 GB RAM per processor.

Figure 6

(a) $T_N$ vs. $U/t$ for system sizes with $4^3,{16}^3$, and ${40}^3$ sites. (b) The spin structure factor $\mathrm{S}(\pi ,\pi ,\pi$) vs. temperature, for system sizes $4^3,6^3,{16}^3$, and ${40}^3$ at $U=8.0t$. The results are obtained using a $4^3$ traveling cluster and 16 processors. The data is averaged over 200 configurations obtained from 2000 MC steps for the cases $N=4^3$ and $N=6^3$, while only 1250 MC steps were used for $N={16}^3$ and $N={40}^3$ with an average over 50 configurations. The magenta arrow indicates the thermodynamic $T_N$ in this case. Panels (c) and (d) present the corresponding results in two dimensions. (c) $T_N$ vs. $U/t$ and (d) the corresponding $\mathrm{S}(\pi ,\pi$)'s for system sizes $8^2,{100}^2$, and ${256}^2$. The results are obtained using a $4^2$ traveling cluster and 16 processors.

Figure 7

Memory and time needed for a single diagonalization using Scalable LAPACK. In (a) we show the RAM (in Gb) per processor that is required to store the arrays necessary to diagonalize a double complex Hermitian matrix. The gray region is defined by a limit of 4Gb RAM per processor, which can diagonalize a Hamiltonian for a ${24}^3$ system (27 648-dimensional matrix) with 8 processors. The data is also presented for 4 and 16 processors; (b) the time needed for a single diagonalization for ${16}^3$ and ${24}^3$ systems vs. the number of processors $N_Q$ used in Scalable LAPACK.

Figure 8

(a) Antiferromagnetic order, $\mathrm{S}(\pi ,\pi$), for the two-dimensional Hubbard model at $U/t=8$ using a $N={32}^2$ lattice with difference traveling cluster sizes as indicated. (b) The one-dimensional checkerboard charge order parameter vs. temperature for the adiabatic Holstein model. Panel (c) shows the corresponding average energy of the system with temperature. The data is presented for two schemes indicated as “central” and “corner”; see text for discussion. The parameters for panels (b) and (c) are the same. The results shown are for $K/t=1$.